"""
核反应堆子通道分析程序 - 基于两流体模型

本程序实现了一个简化的核反应堆子通道分析工具，用于模拟和分析燃料棒束中冷却剂的流动特性。
采用两流体模型处理气液两相流动，考虑了轴向和横向的质量、动量传递。

主要功能：
1. 模拟子通道内的两相流动
2. 计算空泡率、压力分布
3. 分析液相和气相的速度场
4. 考虑子通道间的横向混合

物理模型：
- 两流体模型：分别处理气相和液相
- 考虑相间作用力
- 包含重力和压力梯度效应
- 简化的横向混合模型

数值方法：
- 采用显式时间推进格式
- 交错网格离散
- 迭代求解直至收敛
- 使用稀疏矩阵求解线性方程组
作者：[李泉伯 2024-12-12]
"""

import numpy as np
from CoolProp.CoolProp import PropsSI
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import os
import csv
from matplotlib.animation import FuncAnimation
import multiprocessing as mp
from functools import partial
from scipy.sparse import csr_matrix, lil_matrix
from scipy.sparse.linalg import gmres
import numpy.linalg as la
from numba import jit, njit, prange
import numba as nb
import json
import traceback
import scipy.optimize as opt


# 设置matplotlib的中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']  
plt.rcParams['axes.unicode_minus'] = False  


class SubchannelAnalysis:
    """
    子通道分析的主类，实现了两流体模型的求解
    
    该类实现了基于两流体模型的核反应堆子通道分析，可以模拟和分析3×3棒束16子通道
    系统中的冷却剂流动特性。
    
    主要属性：
        n_axial (int): 轴向节点数
        n_channels (int): 子通道数量
        channel_connections (dict): 子通道连接关系
        connection_coefficients (dict): 子通道间的连接系数
        
    主要方法：
        solve_mass_momentum(): 求解质量和动量守恒方程
        calculate_cross_flow(): 计算横向流动
        get_fluid_properties(): 获取流体物性
        run_simulation(): 运行完整的模拟计算
    """
    
    def __init__(self, n_axial=20):
        """
        初始化3×3棒束16子通道分析器
        
        参数：
            n_axial (int): 轴向网格节点数，默认为2000
            
        初始化的主要内容：
            1. 几何参数设置
            2. 工况参数定义
            3. 网格划分
            4. 场变量初始化
            5. 子通道连接关系建立
        """
        # === 基本几何参数（单位：米）===
        self.rod_diameter = 14.5e-3  # 燃料棒直径
        self.rod_pitch = 18.7e-3    # 燃料棒间距
        self.bundle_size = 3        # 3×3棒束
        self.n_channels = 16        # 总子通道数
        self.n_axial = n_axial      # 轴向网格数
        self.channel_length = 4.0   # 通道总长度 (m)
        self.dz = self.channel_length / n_axial  # 轴向网格间距 (m)
        
        # === 入口工况参数 ===
        self.inlet_pressure = 15.5e6    # 入口压力 (Pa)
        self.inlet_temp = 583.15        # 入口温度 (K)
        self.inlet_mass_flux = 3500.0   # 入口质量流量 (kg/m²·s)
        self.outlet_pressure = 15.3e6   # 出口压力 (Pa)
        self.average_heat_flux = 1.2e6  # 平均热流密度 (W/m²)
        
        # === 时间步长相关参数 ===
        self.dt = 1e-4  # 时间步长 (s)
        self.current_time = 0.0  # 当前时间 (s)
        self.save_interval = 100  # 保存间隔
        
        # === 定义子通道类型（16子通道布局）===
        self.channel_types = {
            1: 'corner',   # 左上角
            2: 'edge',     # 上边
            3: 'edge',     # 上边
            4: 'corner',   # 右上角
            5: 'edge',     # 左边
            6: 'center',   # 中心
            7: 'center',   # 中心
            8: 'edge',     # 右边
            9: 'edge',     # 左边
            10: 'center',  # 中心
            11: 'center',  # 中心
            12: 'edge',    # 右边
            13: 'corner',  # 左下角
            14: 'edge',    # 下边
            15: 'edge',    # 下边
            16: 'corner'   # 右下角
        }
        
        # === 子通道几何参数定义 ===
        self.channel_geometry = {
            'corner': {
                'flow_area': 5.050e-5,        # 流通面积 (m²)
                'wetted_perimeter': 2.835e-2,  # 湿周长 (m)
                'heated_perimeter': 1.14e-2    # 加热周长 (m)
            },
            'edge': {
                'flow_area': 1.177e-4,        # 流通面积 (m²)
                'wetted_perimeter': 4.149e-2,  # 湿周长 (m)
                'heated_perimeter': 2.28e-2    # 加热周长 (m)
            },
            'center': {
                'flow_area': 1.868e-4,        # 流通面积 (m²)
                'wetted_perimeter': 4.548e-2,  # 湿周长 (m)
                'heated_perimeter': 4.56e-2    # 加热周长 (m)
            }
        }
        
        # === 计算总流通面积和质量流量 ===
        self.total_flow_area = sum(self.channel_geometry[self.channel_types[j+1]]['flow_area'] 
                                  for j in range(self.n_channels))
        
        # === 定义子通道连接关系（16子通道拓扑结构）===
        self.channel_connections = {
            1: [2, 5],           # 1号通道连接2,5
            2: [1, 3, 6],        # 2号通道连接1,3,6
            3: [2, 4, 7],        # 3号通道连接2,4,7
            4: [3, 8],           # 4号通道连接3,8
            5: [1, 6, 9],        # 5号通道连接1,6,9
            6: [2, 5, 7, 10],    # 6号通道连接2,5,7,10
            7: [3, 6, 8, 11],    # 7号通道连接3,6,8,11
            8: [4, 7, 12],       # 8号通道连接4,7,12
            9: [5, 10, 13],      # 9号通道连接5,10,13
            10: [6, 9, 11, 14],  # 10号通道连接6,9,11,14
            11: [7, 10, 12, 15], # 11号通道连接7,10,12,15
            12: [8, 11, 16],     # 12号通道连接8,11,16
            13: [9, 14],         # 13号通道连接9,14
            14: [13, 15, 10],    # 14号通道连接13,15,10
            15: [14, 16, 11],    # 15号通道连接14,16,11
            16: [15, 12]         # 16号通道连接15,12
        }
        
        # === 计算连接系数（gap宽度）===
        self.connection_gaps = {}
        for i in range(1, 17):
            for j in self.channel_connections[i]:
                if (i, j) not in self.connection_gaps and (j, i) not in self.connection_gaps:
                    # 根据子通道类型确定连接宽度
                    if self.channel_types[i] == 'corner' or self.channel_types[j] == 'corner':
                        gap_width = 0.02  # 角通道连接宽度 (m)
                    else:
                        gap_width = 0.025 # 其他通道连接宽度 (m)
                    self.connection_gaps[(i, j)] = gap_width
                    self.connection_gaps[(j, i)] = gap_width
        
        # === 创建结果保存文件夹 ===
        self.results_folder = 'transient_results'
        if not os.path.exists(self.results_folder):
            os.makedirs(self.results_folder)
        
        # === 初始化场变量数组 ===
        self.initialize_field_variables()
        
        # === 入口工况参数设置 ===
        self.inlet_pressure = 15.5e6  # 入口压力 (Pa)
        self.inlet_temp = 300         # 入口温度 (K)
        self.mass_flux = 1500         # 质量流量 (kg/m²s)
        
        # === 场变量初始化 ===
        # 空泡率场初始化
        self.void_fraction = np.zeros((n_axial, self.n_channels))    # 空泡率 (-)
        self.void_fraction[0, :] = 0.1  # 入口处初始空泡率设为10%

        # 速度场初始化
        self.liquid_velocity = np.zeros((n_axial, self.n_channels))  # 液相速度 (m/s)
        self.vapor_velocity = np.zeros((n_axial, self.n_channels))   # 气相速度 (m/s)

        # 根据质量流量计算入口速度
        rho_l = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')['density']
        self.liquid_velocity[0, :] = self.mass_flux / rho_l * (1 - self.void_fraction[0, :])  
        self.vapor_velocity[0, :] = self.liquid_velocity[0, :] * 0.1  # 气相速度初始为液相的10%

        # 压力和温度场初始化
        self.pressure = np.full((n_axial, self.n_channels), self.inlet_pressure)  # 压力场 (Pa)
        self.temperature = np.full((n_axial, self.n_channels), self.inlet_temp)   # 温度场 (K)
        self.liquid_enthalpy = np.zeros((n_axial, self.n_channels))  # 液相焓 (J/kg)
        self.vapor_enthalpy = np.zeros((n_axial, self.n_channels))   # 气相焓 (J/kg)
        
        # === 两流体模型物理量初始化 ===
        # 相分数
        self.alpha_l = np.zeros((n_axial, self.n_channels))  # 液相分数 (-)
        self.alpha_v = np.zeros((n_axial, self.n_channels))  # 气相分数 (-)
        
        # 速度场
        self.u_l = np.zeros((n_axial, self.n_channels))  # 液相轴向速度 (m/s)
        self.u_v = np.zeros((n_axial, self.n_channels))  # 气相轴向速度 (m/s)
        self.v_l = np.zeros((n_axial, self.n_channels))  # 液相横向速度 (m/s)
        self.v_v = np.zeros((n_axial, self.n_channels))  # 气相横向速度 (m/s)
        self.w_l = np.zeros((n_axial, self.n_channels))  # 液相z方向速度 (m/s)
        self.w_v = np.zeros((n_axial, self.n_channels))  # 气相z方向速度 (m/s)
        
        # 焓和温度场
        self.h_l = np.zeros((n_axial, self.n_channels))  # 液相焓 (J/kg)
        self.h_v = np.zeros((n_axial, self.n_channels))  # 气相焓 (J/kg)
        self.T_l = np.zeros((n_axial, self.n_channels))  # 液相温度 (K)
        self.T_v = np.zeros((n_axial, self.n_channels))  # 气相温度 (K)
        self.T_c = np.zeros((n_axial, self.n_channels))  # 包壳温度 (K)
        
        # 相变和传热参数
        self.gamma = np.zeros((n_axial, self.n_channels))  # 相变质量流量 (kg/m³·s)
        self.q_wl = np.zeros((n_axial, self.n_channels))  # 液相壁面热流密度 (W/m²)
        self.q_wv = np.zeros((n_axial, self.n_channels))  # 气相壁面热流密度 (W/m²)
        
        # === 瞬态计算参数 ===
        self.total_time = 5.0         # 总计算时间 (s)
        self.time_step = 0.001        # 时间步长 (s)
        self.current_time = 0.0       # 当前时间 (s)
        self.max_iterations = 1000    # 每个时间步内的最大迭代次数
        self.convergence_criteria = 1e-6  # 收敛判据
        
        # === 保存上一时间步的场值 ===
        self.prev_void_fraction = np.zeros((n_axial, self.n_channels))
        self.prev_liquid_velocity = np.zeros((n_axial, self.n_channels))
        self.prev_vapor_velocity = np.zeros((n_axial, self.n_channels))
        self.prev_pressure = np.zeros((n_axial, self.n_channels))
        self.prev_temperature = np.zeros((n_axial, self.n_channels))
        
        # === 初始化物理参数 ===
        self._initialize_physical_parameters()
        
        # === 结果保存参数 ===
        self.save_interval = 100  # 保存间隔步数
        self.results_folder = 'transient_results'  # 结果保存文件夹
        
        # 初始化所有场变量
        self.void_fraction = np.zeros((n_axial, self.n_channels))
        self.pressure = np.full((n_axial, self.n_channels), self.inlet_pressure)
        self.temperature = np.full((n_axial, self.n_channels), self.inlet_temp)
        self.liquid_velocity = np.zeros((n_axial, self.n_channels))
        self.vapor_velocity = np.zeros((n_axial, self.n_channels))
        self.T_l = np.full((n_axial, self.n_channels), self.inlet_temp)
        self.T_v = np.full((n_axial, self.n_channels), self.inlet_temp)
        self.T_c = np.full((n_axial, self.n_channels), self.inlet_temp)
        self.q_wl = np.zeros((n_axial, self.n_channels))
        self.q_wv = np.zeros((n_axial, self.n_channels))
        self.gamma = np.zeros((n_axial, self.n_channels))
        
        # 初始化历史数据字典
        self.history = {
            # 时间相关
            'time': [],
            'time_steps': [],
            
            # 场变量历史
            'void_fraction_history': [],
            'pressure_history': [],
            'temperature_history': [],
            'liquid_velocity_history': [],
            'vapor_velocity_history': [],
            'liquid_temp_history': [],
            'vapor_temp_history': [],
            'wall_temp_history': [],
            'heat_flux_liquid_history': [],
            'heat_flux_vapor_history': [],
            'phase_change_rate_history': [],
            
            # 统计量历史
            'max_void_fraction': [],
            'avg_pressure': [],
            'max_temperature': [],
            'max_liquid_velocity': [],
            'max_vapor_velocity': [],
            'residuals': [],
            'iterations': []
        }
        
        # 在初始化部分添加 w_prime_ij 数组
        self.w_prime_ij = np.zeros((self.n_channels, self.n_channels))  # 湍流混合系数矩阵
        
        # 初始化横向流动和混合系数矩阵
        self.w_ij = np.zeros((self.n_channels, self.n_channels))  # 横向流动系数矩阵
        
        # 计算通道直径
        self.channel_diameter = 2 * np.sqrt(self.channel_geometry[self.channel_types[1]]['flow_area'] / np.pi)
        
        # === 添加界面张力模型参数 ===
        self.surface_tension = 0.072  # 水的表面张力系数 (N/m)
        self.contact_angle = 45  # 接触角 (度)
        
        # === 添加湍流模型参数 (k-ε模型) ===
        self.k = np.zeros((n_axial, self.n_channels))  # 湍流动能
        self.epsilon = np.zeros((n_axial, self.n_channels))  # 湍流耗散率
        self.C_mu = 0.09  # k-ε模型常数
        self.C_1 = 1.44  # k-ε模型常数
        self.C_2 = 1.92  # k-ε模型常数
        self.sigma_k = 1.0  # k的湍流Schmidt数
        self.sigma_epsilon = 1.3  # ε的湍流Schmidt数
        
        # === 添加热耦合效应参数 ===
        self.T_fuel = np.zeros((n_axial, self.n_channels))  # 燃料温度场
        self.T_clad = np.zeros((n_axial, self.n_channels))  # 包壳温度场
        self.k_fuel = 3.0  # 燃料导热系数 (W/m·K)
        self.k_clad = 16.0  # 包壳导热系数 (W/m·K)
        self.gap_conductance = 10000.0  # 间隙导热系数 (W/m²·K)

        # 添加数值计算控制参数
        self.min_dt = 1e-8  # 最小时间步长
        self.max_dt = 1e-2  # 最大时间步长
        self.max_gamma = 1e3  # 最大相变率
        self.dynamic_tolerance = 1e-4  # 动态调整时间步长的误差阈值

    def initialize_field_variables(self):
        """
        初始化16子通道的场变量
        
        该方法初始化所有需要的物理��变量数组，包括:
        - 空泡率场
        - 压力场
        - 温度场
        - 速度场
        - 焓场
        - 传热相关场
        
        所有数组的大小为(n_axial, n_channels)
        """
        # 创建场变量数组
        self.void_fraction = np.zeros((self.n_axial, self.n_channels))
        self.pressure = np.zeros((self.n_axial, self.n_channels))
        self.temperature = np.zeros((self.n_axial, self.n_channels))
        self.liquid_velocity = np.zeros((self.n_axial, self.n_channels))
        self.vapor_velocity = np.zeros((self.n_axial, self.n_channels))
        self.T_l = np.zeros((self.n_axial, self.n_channels))
        self.T_v = np.zeros((self.n_axial, self.n_channels))
        self.T_c = np.zeros((self.n_axial, self.n_channels))
        self.h_l = np.zeros((self.n_axial, self.n_channels))
        self.h_v = np.zeros((self.n_axial, self.n_channels))
        self.q_wl = np.zeros((self.n_axial, self.n_channels))
        self.q_wv = np.zeros((self.n_axial, self.n_channels))
        self.gamma = np.zeros((self.n_axial, self.n_channels))

    def _initialize_physical_parameters(self):
        """
        初始化物理参数，考虑16子通道的特性
        
        该方法完成以下初始化工作�����
        1. 设置入口边界条件
        2. 初始化��向分布
        3. 为��同类型的子通道设置合适的初始值
        4. 考虑重力效应
        
        注意：
            - 角通道使用0.7的流量系数
            - 边通道使用0.9的流量系数
            - 中心通道使用1.0的流量系数
        """
        # 设置入口边界条件
        for j in range(self.n_channels):
            channel_type = self.channel_types[j+1]
            flow_area = self.channel_geometry[channel_type]['flow_area']
            
            # 入口空泡率和相速度
            self.void_fraction[0, j] = 0.0  # 入口无蒸汽
            self.alpha_l[0, j] = 1.0
            self.alpha_v[0, j] = 0.0
            
            # 入口压力和温度
            self.pressure[0, j] = self.inlet_pressure
            self.temperature[0, j] = self.inlet_temp
            self.T_l[0, j] = self.inlet_temp
            self.T_v[0, j] = self.inlet_temp
            self.T_c[0, j] = self.inlet_temp
            
            # 根据子通道类型计算入口速度
            rho_l = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')['density']
            # 根据子通道类型调整质量流量
            if channel_type == 'corner':
                mass_flux_factor = 0.7  # 角通道流量系数
            elif channel_type == 'edge':
                mass_flux_factor = 0.9  # 边通道流量系数
            else:  # center
                mass_flux_factor = 1.0  # 中心通道流量系数
                
            # 计算入口速度
            self.liquid_velocity[0, j] = (self.inlet_mass_flux * mass_flux_factor) / (rho_l * flow_area)
            self.vapor_velocity[0, j] = 0.0
            
            # 计算入口焓
            props_l = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')
            props_v = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'vapor')
            self.h_l[0, j] = props_l['enthalpy']
            self.h_v[0, j] = props_v['enthalpy']

        # 初始化轴向分布
        for i in range(1, self.n_axial):
            # 使用入口条件初始化整个流道
            self.void_fraction[i, :] = self.void_fraction[0, :]
            self.pressure[i, :] = self.pressure[0, :] - i * self.dz * 9.81 * 1000  # 考虑重力压降
            self.temperature[i, :] = self.temperature[0, :]
            self.liquid_velocity[i, :] = self.liquid_velocity[0, :]
            self.vapor_velocity[i, :] = self.vapor_velocity[0, :]
            self.T_l[i, :] = self.T_l[0, :]
            self.T_v[i, :] = self.T_v[0, :]
            self.T_c[i, :] = self.T_c[0, :]
            self.h_l[i, :] = self.h_l[0, :]
            self.h_v[i, :] = self.h_v[0, :]

    def calculate_initial_mass_flux_distribution(self):
        """
        计算初始质量流量分布
        
        根据子通道类型计算每个子通道的初始质量流量分布。
        考虑不同类型子通道的流动特性，使用不同的流量系数。
        
        返回：
            numpy.ndarray: 每个子通道的质量流量分布 (kg/m²·s)
        """
        total_flow_area = 0
        for j in range(self.n_channels):
            channel_type = self.channel_types[j+1]
            total_flow_area += self.channel_geometry[channel_type]['flow_area']
        
        # 计算平均质量流量
        avg_mass_flux = self.total_mass_flow / total_flow_area
        
        # 为每个子通道分配质量流量
        mass_flux_distribution = np.zeros(self.n_channels)
        for j in range(self.n_channels):
            channel_type = self.channel_types[j+1]
            if channel_type == 'corner':
                mass_flux_distribution[j] = avg_mass_flux * 0.7
            elif channel_type == 'edge':
                mass_flux_distribution[j] = avg_mass_flux * 0.9
            else:  # center
                mass_flux_distribution[j] = avg_mass_flux * 1.0
        
        return mass_flux_distribution

    def initialize_heat_flux(self):
        """
        初始化各子通道热流密度分布
        
        该方法：
        1. 考虑轴向功率分布（采用余弦分布）
        2. 考虑径向功率分布（不同类型子通道使用不同的功率因子）
        3. 初始时假设所有热量都传给液相
        
        功率因子：
        - 角通道：0.8
        - 边通道：0.9
        - 中心通道：1.0
        """
        for j in range(self.n_channels):
            channel_type = self.channel_types[j+1]
            heated_perimeter = self.channel_geometry[channel_type]['heated_perimeter']
            
            # 考虑轴向功率分布
            for i in range(self.n_axial):
                z = i * self.dz
                # 余弦分布的轴向功率因子
                axial_power_factor = np.cos(np.pi * z / (self.n_axial * self.dz))
                
                # 根据子通道类型设置径向功率因子
                if channel_type == 'corner':
                    radial_power_factor = 0.8
                elif channel_type == 'edge':
                    radial_power_factor = 0.9
                else:  # center
                    radial_power_factor = 1.0
                    
                # 计算局部热流密度
                q_base = self.average_heat_flux * axial_power_factor * radial_power_factor
                self.q_wl[i, j] = q_base * heated_perimeter
                self.q_wv[i, j] = 0.0  # 初始时假设所有热量都传给液相

    def get_fluid_properties(self, p, t, phase):
        """获取指定状态点的流体物性"""
        try:
            # 确保参数在CoolProp的有效范围内
            p = max(611.655, min(p, 2.2064e7))  # 水的三相点到临界点压力范围
            t = max(273.16, min(t, 647.096))    # 水的三相点到临界点温度范围
            
            if phase == 'liquid':
                return {
                    'density': PropsSI('D', 'P', p, 'T', t, 'Water'),
                    'viscosity': PropsSI('V', 'P', p, 'T', t, 'Water'),
                    'enthalpy': PropsSI('H', 'P', p, 'T', t, 'Water'),
                    'specific_heat': PropsSI('C', 'P', p, 'T', t, 'Water'),
                    'thermal_conductivity': PropsSI('L', 'P', p, 'T', t, 'Water')
                }
            else:  # vapor
                return {
                    'density': PropsSI('D', 'P', p, 'Q', 1, 'Water'),
                    'viscosity': PropsSI('V', 'P', p, 'Q', 1, 'Water'),
                    'enthalpy': PropsSI('H', 'P', p, 'Q', 1, 'Water'),
                    'specific_heat': PropsSI('C', 'P', p, 'Q', 1, 'Water'),
                    'thermal_conductivity': PropsSI('L', 'P', p, 'Q', 1, 'Water')
                }
        except Exception as e:
            print(f"警告：物性计算失败 (P={p}, T={t})")
            print(f"错误信息: {str(e)}")
            return None
            
    def calculate_cross_flow(self, i, j, k):
        """
        计算子通道j和k之间的横向流动
        
        基于压力差和几何特征计算相邻子通道间的横向质量流量。
        考虑不同类型子通道的阻力特性。
        
        参数:
            i (int): 轴向位置索引
            j (int): 源子通道编号
            k (int): 目标子通道编号
            
        返回:
            float: 横向质量流量 (kg/s)，正值表示从j到k的流动
            
        注意:
            - 角通道的阻力系数为1.5
            - 边通道的阻力系数为1.2
            - 中心通道的阻力系数为1.0
        """
        try:
            # 检查子通道连接是否存在
            if k not in self.channel_connections[j]:
                return 0.0
            
            # 获取连接宽度
            gap_width = self.connection_gaps.get((j, k), 0.0)
            
            # 计算压力差
            dp = self.pressure[i, j-1] - self.pressure[i, k-1]
            
            # 获取子通道类型
            j_type = self.channel_types[j]
            k_type = self.channel_types[k]
            
            # 计算横向流动阻力系数
            if j_type == 'corner' or k_type == 'corner':
                resistance_factor = 1.5  # 角通道的阻力系数较大
            elif j_type == 'edge' or k_type == 'edge':
                resistance_factor = 1.2  # 边通道的阻力系数中等
            else:
                resistance_factor = 1.0  # 中心通道的阻力系数标准值
            
            # 计算横向流动系数
            s_jk = gap_width / (resistance_factor * self.dz)
            
            # 计算横向质量流量
            w_jk = s_jk * dp
            
            return w_jk
            
        except Exception as e:
            print(f"计算横向流动时出错: {str(e)}")
            return 0.0

    def calculate_mixing_coefficient(self, i, j, k):
        """
        计算子通道j和k之间的湍流混合系数
        
        基于流动特性和子通道类型计算湍流混合系数。
        考虑不同类型子通道的混合特性。
        
        参数：
            i (int): 轴向位置索引
            j (int): 第一个子通道编号
            k (int): 第二个子通道编号
            
        返回：
            float: 湍流混合系数 (无量纲)
            
        注意：
            - 角通道混合系数为基准值的0.7倍
            - 边通道混合系数为基准值的0.85倍
            - 中心通道混合系数为基准值
        """
        try:
            # 获取子通道类型
            j_type = self.channel_types[j]
            k_type = self.channel_types[k]
            
            # 获取平均流速
            v_j = self.liquid_velocity[i, j-1]
            v_k = self.liquid_velocity[i, k-1]
            v_avg = 0.5 * (v_j + v_k)
            
            # 基础混合系数
            beta_0 = 0.02
            
            # 根据子通道类型调整混合系数
            if j_type == 'corner' or k_type == 'corner':
                type_factor = 0.7  # 角通道混合较弱
            elif j_type == 'edge' or k_type == 'edge':
                type_factor = 0.85  # 边通道混合中等
            else:
                type_factor = 1.0  # 中心通道混合正常
            
            # 计算实际混合系数
            beta = beta_0 * type_factor * v_avg
            
            return beta
            
        except Exception as e:
            print(f"计算混合系数时出错: {str(e)}")
            return 0.0

    def calculate_turbulent_mixing(self, i, j, k, phi):
        """
        计算子通道j和k之间的湍流混合传输
        
        基于浓度梯度和混合系数计算相邻子通道间的湍流混合传输。
        
        参数:
            i (int): 轴向位置索引
            j, k (int): 相邻子通道编号
            phi (numpy.ndarray): 被混合的物理量（温度、焓等）
            
        返回:
            float: 混合传输量，正值表示从k到j的传输
            
        注意：
            - 混合方向由浓度梯度决定
            - 混合强度与连接宽度和浓度差值成正比
        """
        try:
            # 检查子通道连接
            if k not in self.channel_connections[j]:
                return 0.0
                
            # 获取连接宽度
            gap_width = self.connection_gaps.get((j, k), 0.0)
            
            # 计算混合系数
            beta = self.calculate_mixing_coefficient(i, j, k)
            
            # 计算物理量差值
            dphi = phi[i, j-1] - phi[i, k-1]
            
            # 计算混合流量
            w_prime = beta * gap_width * abs(dphi)
            
            # 确定混合方向
            if dphi > 0:
                return -w_prime  # 从j到k的混合
            else:
                return w_prime   # 从k到j的混合
                
        except Exception as e:
            print(f"计算湍流混合时出错: {str(e)}")
            return 0.0

    def update_cross_flow_terms(self, i):
        """
        更新所有子通道间的横向流动项
        """
        try:
            # 初始化横向流动数组
            self.w_ij = np.zeros((self.n_channels, self.n_channels))
            self.w_prime_ij = np.zeros((self.n_channels, self.n_channels))
            
            # 遍历所有子通道对
            for j in range(1, self.n_channels + 1):
                for k in self.channel_connections[j]:
                    if k > j:  # 避免重复计算
                        # 计算横向流动
                        w_jk = self.calculate_cross_flow(i, j, k)
                        self.w_ij[j-1, k-1] = w_jk
                        self.w_ij[k-1, j-1] = -w_jk
                        
                        # 计算湍流混合
                        # 温度混合
                        w_prime_t = self.calculate_turbulent_mixing(i, j, k, self.temperature)
                        # 焓混合
                        w_prime_h = self.calculate_turbulent_mixing(i, j, k, self.h_l)
                        
                        # 取平均值作为最终的湍流混合项
                        w_prime = 0.5 * (w_prime_t + w_prime_h)
                        self.w_prime_ij[j-1, k-1] = w_prime
                        self.w_prime_ij[k-1, j-1] = -w_prime
                        
        except Exception as e:
            print(f"更新横向流动项时出错: {str(e)}")

    def solve_mass_momentum(self):
        """求解质量和动量方程"""
        print("\n=== 开始求解质量和动量方程 ===")
        dt = self.dt  # 使用类的时间步长
        max_iterations = self.max_iterations
        tolerance = self.convergence_criteria
        
        for n in range(max_iterations):
            max_error = 0.0
            
            # 保存旧值用于计算误差
            old_void = np.copy(self.void_fraction)
            old_vl = np.copy(self.liquid_velocity)
            old_vv = np.copy(self.vapor_velocity)
            
            print(f"\n迭代步数: {n}")
            print(f"当前最大空泡率: {np.max(self.void_fraction):.4f}")
            print(f"当前最大液相速度: {np.max(self.liquid_velocity):.4f} m/s")
            print(f"当前最大气相速度: {np.max(self.vapor_velocity):.4f} m/s")
            
            # ���新横向流动项
            self.update_cross_flow_terms(n)
            
            for i in range(1, self.n_axial):
                for j in range(self.n_channels):
                    # 获���子通道类型和几何参数
                    channel_type = self.channel_types[j+1]
                    flow_area = self.channel_geometry[channel_type]['flow_area']
                    wetted_perimeter = self.channel_geometry[channel_type]['wetted_perimeter']
                    
                    print(f"\n处理子通道 {j+1} (类型: {channel_type})")
                    print(f"轴向位置: {i*self.dz:.3f} m")
                    print(f"流通面积: {flow_area:.6f} m²")
                    print(f"湿周长: {wetted_perimeter:.6f} m")
                    
                    # 获取流体属性
                    props_l = self.get_fluid_properties(self.pressure[i,j], self.temperature[i,j], 'liquid')
                    props_v = self.get_fluid_properties(self.pressure[i,j], self.temperature[i,j], 'vapor')
                    
                    if props_l is None or props_v is None:
                        print("警告: 无法获取流体物性")
                        continue
                        
                    rho_l = props_l['density']
                    rho_v = props_v['density']
                    
                    print(f"液相密度: {rho_l:.2f} kg/m³")
                    print(f"气相密度: {rho_v:.2f} kg/m³")
                    
                    # 计算横向流动的净贡献
                    net_cross_flow = 0.0
                    net_mixing = 0.0
                    for k in self.channel_connections[j+1]:
                        net_cross_flow += self.w_ij[j, k-1]
                        net_mixing += self.w_prime_ij[j, k-1]
                    
                    print(f"横向流动净贡献: {net_cross_flow:.4f} kg/s")
                    print(f"湍流混合净贡献: {net_mixing:.4f} kg/s")
                    
                    # 计算摩擦压降
                    f_l = self.calculate_friction_factor(self.liquid_velocity[i,j], 
                                                            props_l['viscosity'],
                                                        flow_area,
                                                        wetted_perimeter)
                                                    
                    f_v = self.calculate_friction_factor(self.vapor_velocity[i,j],
                                                            props_v['viscosity'],
                                                        flow_area,
                                                        wetted_perimeter)
                    
                    print(f"液相摩擦因子: {f_l:.4f}")
                    print(f"气相摩擦因子: {f_v:.4f}")
                    
                    # 使用原有的数值方法更新场变量
                    alpha = self.void_fraction[i-1,j]
                    vl = self.liquid_velocity[i-1,j]
                    vv = self.vapor_velocity[i-1,j]
                    
                    # 空泡率方程
                    new_alpha = alpha + dt * (
                        -vv * (alpha - self.void_fraction[i-1,j]) / self.dz
                        + net_cross_flow / (rho_v * flow_area)
                    )
                    new_alpha = np.clip(new_alpha, 0, 0.9999)
                    
                    print(f"更新后的空泡率: {new_alpha:.4f}")
                    
                    # 液相动量方程
                    new_vl = vl + dt * (
                        -vl * (vl - self.liquid_velocity[i-1,j]) / self.dz
                        - (1/rho_l) * (self.pressure[i,j] - self.pressure[i-1,j]) / self.dz
                        - f_l * vl * abs(vl) / (2 * self.channel_diameter)
                        - 9.81
                    )
                    
                    print(f"更新后的液相速度: {new_vl:.4f} m/s")
                    
                    # 气相动量方程
                    new_vv = vv + dt * (
                        -vv * (vv - self.vapor_velocity[i-1,j]) / self.dz
                        - (1/rho_v) * (self.pressure[i,j] - self.pressure[i-1,j]) / self.dz
                        - f_v * vv * abs(vv) / (2 * self.channel_diameter)
                        - 9.81
                    )
                    
                    print(f"更新后的气相速度: {new_vv:.4f} m/s")
                    
                    # 更新场变量
                    self.void_fraction[i,j] = new_alpha
                    self.liquid_velocity[i,j] = new_vl
                    self.vapor_velocity[i,j] = new_vv
                    
                    # 计算最大误差
                    error = max(
                        abs(new_alpha - old_void[i,j]),
                        abs(new_vl - old_vl[i,j]),
                        abs(new_vv - old_vv[i,j]))
                    max_error = max(max_error, error)
            
            print(f"\n当前迭代最大误差: {max_error:.6f}")
            
            # 检查收敛性并动态调整时间步长
            if max_error > self.dynamic_tolerance:
                dt *= 0.5
                print(f"误差过大，时间步长调整为 {dt:.2e}")
                if dt < self.min_dt:
                    print("错误：时间步长过小，计算终止")
                    break
            elif max_error < tolerance:
                print(f"\n计算在第{n}次迭代后收敛")
                break
            else:
                dt = min(dt * 1.1, self.max_dt)  # 适当增加时间步长但不超过最大值
            
            self.dt = dt
    
    def run_simulation(self):
        """运行完整的模拟计算"""
        print("开始数值计算...")
        
        time = 0
        step = 0
        
        while time < self.total_time:
            print(f"当前时间：{time:.3f} s")
            
            # 使用自适应时间步长解
            success = self.solve_two_fluid_equations()
            
            if not success:
                print(f"警告：在时间 {time:.3f} s 处计算不稳定")
                self.time_step *= 0.1
                if self.time_step < 1e-8:
                    print("错误：时间步长过小，计算终止")
                    break
                continue
            
            # 更新时间
            time += self.time_step
            step += 1
            
            # 保存结果
            if step % self.save_interval == 0:
                self.save_timestep_results(step)
                
        print("计算完成")
    
    def plot_results(self):
        """绘制计算结果的可视化图，增加错误检查和空值处理"""
        try:
            # 创建网格点
            x = np.arange(self.n_axial)
            y = np.arange(self.n_channels)
            X, Y = np.meshgrid(x, y)
            
            # 检查数据是否为空
            if np.all(self.void_fraction == 0) or np.all(self.pressure == 0):
                print("警告：检测到空数据，请确保计算已完成")
                return
                
            # 创建两个图形窗口
            fig1 = plt.figure(figsize=(20, 15))
            fig2 = plt.figure(figsize=(20, 15))
            
            # 第1个图形窗口：液相参数
            axes1 = []
            for i in range(6):
                ax = fig1.add_subplot(2, 3, i+1, projection='3d')
                axes1.append(ax)
                
            # 第二个图形窗口：气相参数
            axes2 = []
            for i in range(6):
                ax = fig2.add_subplot(2, 3, i+1, projection='3d')
                axes2.append(ax)
                
            # 绘制液相参数
            plots = [
                (self.alpha_l, '液相分数分布'),
                (self.u_l, '液相轴向速度分布 (m/s)'),
                (self.v_l, '液相横向速度分布 (m/s)'),
                (self.T_l, '液相温度分布 (K)'),
                (self.pressure, '压力分布 (Pa)'),
                (self.q_wl if hasattr(self, 'q_wl') else np.zeros_like(self.pressure), 
                 '液相壁面热流密度 (W/m²)')
            ]
            
            for ax, (data, title) in zip(axes1, plots):
                surf = ax.plot_surface(X, Y, data.T, cmap=cm.viridis)
                ax.set_title(title)
                ax.set_xlabel('轴向位置')
                ax.set_ylabel('子通道编号')
                fig1.colorbar(surf, ax=ax, shrink=0.5, aspect=5)
                
            # 绘制气相参数
            plots = [
                (self.alpha_v, '气相分数分布'),
                (self.u_v, '气相轴向速度分布 (m/s)'),
                (self.v_v, '气相横向速度分布 (m/s)'),
                (self.T_v, '气相温度分布 (K)'),
                (self.gamma if hasattr(self, 'gamma') else np.zeros_like(self.pressure), 
                 '相变率 (kg/m³·s)'),
                (self.q_wv if hasattr(self, 'q_wv') else np.zeros_like(self.pressure), 
                 '气相壁面热流密度 (W/m²)')
            ]
            
            for ax, (data, title) in zip(axes2, plots):
                surf = ax.plot_surface(X, Y, data.T, cmap=cm.viridis)
                ax.set_title(title)
                ax.set_xlabel('轴向位置')
                ax.set_ylabel('子通道编号')
                fig2.colorbar(surf, ax=ax, shrink=0.5, aspect=5)
                
            # 设置视角和网格
            for ax in axes1 + axes2:
                ax.view_init(elev=20, azim=45)
                ax.grid(True)
                
            plt.tight_layout()
            plt.show()
            
        except Exception as e:
            print(f"绘图时出错: {str(e)}")

    def solve_two_fluid_equations(self):
        """求解双流体模型方程组 - 添加自适应时间步长控制"""
        # 准备输入数据
        dx = dy = self.dz  # 假设网格等间距
        max_attempts = 5  # 最大重试次数
        
        for attempt in range(max_attempts):
            # 调用Numba优化的求解器
            results = solve_two_fluid_equations_numba(
                self.n_axial, self.n_channels,
                self.alpha_l, self.alpha_v,
                self.rho_l, self.rho_v,
                self.u_l, self.u_v,
                self.v_l, self.v_v,
                self.w_l, self.w_v,
                self.pressure, self.T_l, self.T_v, self.T_c,
                self.h_l, self.h_v,
                dx, dy, self.dz, self.time_step
            )
            
            # 解包结果
            if len(results) == 10:  # 正常返回
                (new_alpha_l, new_alpha_v, new_u_l, new_u_v,
                 new_v_l, new_v_v, new_h_l, new_h_v,
                 new_dt, needs_restart) = results
                
                if needs_restart:
                    print(f"时间步长调整为 {new_dt:.2e}，重新计算...")
                    self.time_step = new_dt
                    continue
                    
                # 更新时间步长
                self.time_step = new_dt
                
                # 更新场值
                self.alpha_l = new_alpha_l
                self.alpha_v = new_alpha_v
                self.u_l = new_u_l
                self.u_v = new_u_v
                self.v_l = new_v_l
                self.v_v = new_v_v
                self.h_l = new_h_l
                self.h_v = new_h_v
                
                print(f"计算成功，当前时间步长: {self.time_step:.2e}")
                return True
                
        print("警告：达到最大重试次数，计算可能不稳定")
        return False

    def calculate_phase_change(self, i, j):
        """计算相变质量流量"""
        try:
            T_c = self.T_c[i,j]  # 包壳温度
            T_sat = self.get_saturation_temperature(self.pressure[i,j])
            T_l = self.T_l[i,j]  # 液相温度
            
            # 过冷沸腾
            if T_l < T_sat and T_c > T_sat:
                h_l = self.calculate_subcooled_boiling_htc(T_c, T_sat, T_l, self.pressure[i,j])
                Q_l = np.pi * self.channel_diameter * h_l * (T_c - T_l)
                gamma_l_sc = Q_l / self.get_latent_heat(self.pressure[i,j])
                gamma_l_sc = np.clip(gamma_l_sc, -self.max_gamma, 0)  # 限制相变率范围
                return gamma_l_sc
                
            # 饱和核态沸腾
            elif T_l >= T_sat and T_c > T_sat:
                h_l = self.calculate_nucleate_boiling_htc(T_c, T_sat, self.pressure[i,j])
                Q_l = np.pi * self.channel_diameter * h_l * (T_c - T_sat)
                gamma_l_sh = Q_l / self.get_latent_heat(self.pressure[i,j])
                gamma_l_sh = np.clip(gamma_l_sh, 0, self.max_gamma)  # 限制相变率范围
                return gamma_l_sh
                
            return 0.0
            
        except Exception as e:
            print(f"计算相变率时出错: {str(e)}")
            return 0.0

    def calculate_turbulent_mixing(self, i, j, phase):
        """计算湍流混合和空泡漂移质量流量"""
        C_e = 0.01  # 涡流扩散率
        L_tur = self.channel_diameter  # 湍流混合长度
        
        mixing_mass = 0.0
        for k in self.channel_connections[j+1]:
            L_gap = self.get_gap_width(j+1, k)  # 获取间隙宽度
            
            if phase == 'liquid':
                rho_diff = self.get_fluid_properties(self.pressure[i,j], 
                    self.temperature[i,j], 'liquid')['density']
                alpha_diff = self.alpha_l[i,j] - self.alpha_l[i,k-1]
                alpha_eq_diff = self.calculate_equilibrium_void_fraction(i,j) - \
                               self.calculate_equilibrium_void_fraction(i,k-1)
            else:
                rho_diff = self.get_fluid_properties(self.pressure[i,j],
                    self.temperature[i,j], 'vapor')['density']
                alpha_diff = self.alpha_v[i,j] - self.alpha_v[i,k-1]
                alpha_eq_diff = (1-self.calculate_equilibrium_void_fraction(i,j)) - \
                               (1-self.calculate_equilibrium_void_fraction(i,k-1))
                
            mixing_mass += C_e * L_gap * L_tur * rho_diff * (alpha_diff - alpha_eq_diff)
            
        return mixing_mass

    def calculate_wall_friction(self, i, j, phase, direction='axial'):
        """计算壁面摩擦力
        
        参数:
        i, j: 网格位置
        phase: 'liquid'或'vapor'
        direction: 'axial'或'lateral'
        """
        if phase == 'liquid':
            rho = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')['density']
            mu = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')['viscosity']
            if direction == 'axial':
                v = self.u_l[i,j]
            else:
                v = self.v_l[i,j]
        else:
            rho = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')['density']
            mu = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')['viscosity']
            if direction == 'axial':
                v = self.u_v[i,j]
            else:
                v = self.v_v[i,j]
            
        # 计算Reynolds数
        Re = rho * abs(v) * self.channel_diameter / mu
        
        # 计算摩擦系数
        if Re < 2300:  # 层流
            f = 64 / Re
        else:  # 湍流
            f = 0.316 * Re**(-0.25)
            
        # 计算壁面摩擦力
        tau_w = 0.5 * f * rho * v * abs(v)
        
        # 考虑两相流修正系数
        if phase == 'liquid':
            phi_l = (1 - self.alpha_l[i,j])**(-1.75)  # 两相流修正系数
            tau_w *= phi_l
            
        return tau_w * np.pi * self.channel_diameter / self.channel_area

    def calculate_interfacial_force(self, i, j, direction='axial'):
        """计算界面剪切力 (eq.19)"""
        rho_l = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')['density']
        
        if direction == 'axial':
            v_rel = self.u_v[i,j] - self.u_l[i,j]  # 相对速度
        else:
            v_rel = self.v_v[i,j] - self.v_l[i,j]
            
        # 计算界面面积密度
        a_i = 6 * self.alpha_v[i,j] / self.channel_diameter  # 假设球形气泡
        
        # 计算界面阻力系数
        C_i = 0.44  # 球形气泡的阻力系数
        
        # 计算界面剪切力
        tau_i = 0.5 * C_i * rho_l * v_rel * abs(v_rel) * a_i
        
        return tau_i

    def calculate_wall_heat_flux(self, i, j, phase):
        """计算壁面热流密度 (eq.6)"""
        if phase == 'liquid':
            # 根据沸腾状态确定换热系数
            T_sat = self.get_saturation_temperature(self.pressure[i,j])
            
            if self.T_l[i,j] < T_sat and self.T_c[i,j] > T_sat:
                # 过冷核态沸腾
                h = self.calculate_subcooled_boiling_htc(
                    self.T_c[i,j], T_sat, self.T_l[i,j], self.pressure[i,j])
                
                
            elif self.T_l[i,j] >= T_sat and self.T_c[i,j] > T_sat:
                # 饱和核态沸腾
                h = self.calculate_nucleate_boiling_htc(
                    self.T_c[i,j], T_sat, self.pressure[i,j])
                
                
            else:
                # 单相对流换热
                h = self.calculate_single_phase_htc(i, j, 'liquid')
                
            q = h * (self.T_c[i,j] - self.T_l[i,j])
        else:
            # 气相单相对流换热
            h = self.calculate_single_phase_htc(i, j, 'vapor')
            q = h * (self.T_c[i,j] - self.T_v[i,j])
            
        return q * np.pi * self.channel_diameter

    def calculate_subcooled_boiling_htc(self, T_w, T_sat, T_l, P):
        """计算过冷核态沸腾换热系数 (eq.8)"""
        F_gam = 0.5  # 气泡冷凝系数
        h = (1 - F_gam) * np.exp(2*P/1260 * (T_w-T_sat)/(2722*(T_w-T_l)))
        return h

    def calculate_nucleate_boiling_htc(self, T_w, T_sat, P):
        """计算饱和核态沸腾换热系数 (eq.9)"""
        h = np.exp(2*P/1260 * (T_w-T_sat)/2722)
        return h

    def calculate_radial_heat_conduction(self, i, j):
        """计算径向热传导项"""
        # 简化处理,假设均匀网格
        dr = self.channel_diameter / 10  # 径向网格间距
        r = np.linspace(dr/2, self.channel_diameter/2, 10)
        
        # 计算径向温度梯度
        d2T_dr2 = np.zeros_like(r)
        for k in range(1, len(r)-1):
            d2T_dr2[k] = (1/r[k]) * (
                (r[k+1]*self.T_c[i,j] - r[k]*self.T_c[i,j])/dr -
                (r[k]*self.T_c[i,j] - r[k-1]*self.T_c[i,j])/dr
            )  # 添加闭合括号
        
        return np.mean(d2T_dr2)  # 返回平均值

    def calculate_fuel_heat_generation(self, i, j):
        """计算燃料热生成"""
        # 简化处理,假设轴向余弦分布
        q_max = 1e6  # 最大线功率密度 (W/m)
        z = i * self.dz
        H = self.channel_length
        q = q_max * np.cos(np.pi * z / H)
        return q / (np.pi * self.channel_diameter)

    def run_transient_simulation(self):
        """运行瞬态计算"""
        print("开始瞬态计算...")
        
        n_time_steps = int(self.total_time / self.time_step)
        
        for time_step in range(n_time_steps):
            # ���存上一时间步的值
            self.save_previous_step_values()
            
            # 在当前时间步内进行迭代求解
            converged, n_iter, final_residual = self.solve_current_time_step()
            
            # 记录残差和迭代次数
            self.history['residuals'].append(final_residual)
            self.history['iterations'].append(n_iter)
            
            # 保存结果
            self.save_timestep_results(time_step)
            
            if not converged:
                print(f"警告：时间步 {time_step + 1} 未达到收敛")
                
            # 更新当前时间
            self.current_time += self.time_step
            
            # 打印计算进度
            if (time_step + 1) % 100 == 0:
                print(f"计算进度: {(time_step + 1)/n_time_steps*100:.1f}% "
                      f"(当前时间: {self.current_time:.3f} s)")
        
        print("瞬态计算完成")
        print(f"结果已保存在 '{self.results_folder}' 文件夹中")
        
    def save_previous_step_values(self):
        """保存上一时间步的场值"""
        self.prev_void_fraction = np.copy(self.void_fraction)
        self.prev_liquid_velocity = np.copy(self.liquid_velocity)
        self.prev_vapor_velocity = np.copy(self.vapor_velocity)
        self.prev_pressure = np.copy(self.pressure)
        self.prev_temperature = np.copy(self.temperature)

    def solve_current_time_step(self):
        """在当前时间步内进行迭代求解"""
        residuals = []
        
        for iteration in range(self.max_iterations):
            # 保存当前迭代步的值用于收敛性检查
            old_void = np.copy(self.void_fraction)
            old_vl = np.copy(self.liquid_velocity)
            old_vv = np.copy(self.vapor_velocity)
            old_p = np.copy(self.pressure)
            old_t = np.copy(self.temperature)
            
            # 求解质量和动量方程
            self.solve_mass_momentum_equations()
            
            # 求解能量方程
            self.solve_energy_equations()
            
            # 检查收敛性
            max_error = self.check_convergence(
                old_void, old_vl, old_vv, old_p, old_t)
            
            residuals.append(max_error)
            
            if max_error < self.convergence_criteria:
                return True, iteration + 1, max_error
                
        return False, self.max_iterations, residuals[-1]
    
    def save_timestep_results(self, time_step):
        """保存当前时间步的结果"""
        try:
            if time_step % self.save_interval == 0:
                # 1. 更新时间相关数据
                self.history['time'].append(self.current_time)
                self.history['time_steps'].append(time_step)
                
                # 2. 更新场变量历史
                self.history['void_fraction_history'].append(np.copy(self.void_fraction))
                self.history['pressure_history'].append(np.copy(self.pressure))
                self.history['temperature_history'].append(np.copy(self.temperature))
                self.history['liquid_velocity_history'].append(np.copy(self.liquid_velocity))
                self.history['vapor_velocity_history'].append(np.copy(self.vapor_velocity))
                self.history['liquid_temp_history'].append(np.copy(self.T_l))
                self.history['vapor_temp_history'].append(np.copy(self.T_v))
                self.history['wall_temp_history'].append(np.copy(self.T_c))
                self.history['heat_flux_liquid_history'].append(np.copy(self.q_wl))
                self.history['heat_flux_vapor_history'].append(np.copy(self.q_wv))
                self.history['phase_change_rate_history'].append(np.copy(self.gamma))
                
                # 3. 更新统计量历史
                self.history['max_void_fraction'].append(np.max(self.void_fraction))
                self.history['avg_pressure'].append(np.mean(self.pressure))
                self.history['max_temperature'].append(np.max(self.temperature))
                self.history['max_liquid_velocity'].append(np.max(self.liquid_velocity))
                self.history['max_vapor_velocity'].append(np.max(self.vapor_velocity))
                
                # 4. 保存CSV文件
                self.save_to_csv(time_step)
                
                # 5. 定期保存完整状态（每100个时间步）
                if time_step % (self.save_interval * 100) == 0:
                    self.save_full_state(time_step)
                    
        except Exception as e:
            print(f"保存时间步 {time_step} 的结果时出错: {str(e)}")
            print(f"错误类型: {type(e).__name__}")
            print(f"错误位置: {traceback.format_exc()}")
            
    def save_to_csv(self, time_step):
        """保存当前时间步的数据到CSV文件"""
        try:
            filename = os.path.join(self.results_folder, f"timestep_{time_step:06d}.csv")
            
            # 使用UTF-8编码保存
            with open(filename, 'w', newline='', encoding='utf-8') as csvfile:
                writer = csv.writer(csvfile)
                
                # 1. 写入时间步信息
                writer.writerow(['时间步信息'])
                writer.writerow(['当前时间 (s)', f"{self.current_time:.6f}"])
                writer.writerow(['时间步编号', time_step])
                writer.writerow([])
                
                # 2. 写入全局参数
                writer.writerow(['全局参数'])
                writer.writerow(['平均压力 (Pa)', f"{np.mean(self.pressure):.2f}"])
                writer.writerow(['最大温度 (K)', f"{np.max(self.temperature):.2f}"])
                writer.writerow(['最大空泡率', f"{np.max(self.void_fraction):.4f}"])
                writer.writerow([])
                
                # 3. 写入详细数据
                headers = [
                    '轴向位置 (m)',
                    '子通道编号',
                    '空泡率',
                    '压力 (Pa)',
                    '液相温度 (K)',
                    '气相温度 (K)',
                    '壁面温度 (K)',
                    '液相速度 (m/s)',
                    '气相速度 (m/s)',
                    '液相壁面热流密度 (W/m²)',
                    '气相壁面热流密度 (W/m²)',
                    '相变率 (kg/m³·s)'
                ]
                writer.writerow(headers)
                
                # 写入每个网格点的数据
                for i in range(self.n_axial):
                    for j in range(self.n_channels):
                        row_data = [
                            f"{i * self.dz:.6f}",
                            j + 1,
                            f"{self.void_fraction[i, j]:.6f}",
                            f"{self.pressure[i, j]:.2f}",
                            f"{self.T_l[i, j]:.2f}",
                            f"{self.T_v[i, j]:.2f}",
                            f"{self.T_c[i, j]:.2f}",
                            f"{self.liquid_velocity[i, j]:.4f}",
                            f"{self.vapor_velocity[i, j]:.4f}",
                            f"{self.q_wl[i, j]:.2f}",
                            f"{self.q_wv[i, j]:.2f}",
                            f"{self.gamma[i, j]:.6f}"
                        ]
                        writer.writerow(row_data)
                        
            print(f"结果已保存到: {filename}")
            
        except Exception as e:
            print(f"保存CSV文件时出错: {str(e)}")
            
    def save_full_state(self, time_step):
        """保存完整的计算状态"""
        try:
            filename = os.path.join(self.results_folder, f"full_state_{time_step:06d}.npz")
            np.savez_compressed(
                filename,
                time=self.current_time,
                time_step=time_step,
                void_fraction=self.void_fraction,
                pressure=self.pressure,
                temperature=self.temperature,
                liquid_velocity=self.liquid_velocity,
                vapor_velocity=self.vapor_velocity,
                T_l=self.T_l,
                T_v=self.T_v,
                T_c=self.T_c,
                q_wl=self.q_wl,
                q_wv=self.q_wv,
                gamma=self.gamma
            )
            print(f"完整状态已保存到: {filename}")
            
        except Exception as e:
            print(f"保存完整状态时出错: {str(e)}")

    def plot_transient_results(self):
        """绘制瞬态计算结果"""
        try:
            # 创建一个包含4个子图的图形
            fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12)) # 移除多余的括号
            
            # 1. 残差历史
            ax1.semilogy(self.history['residuals'])
            ax1.set_title('残差历史')
            ax1.set_xlabel('时间步')
            ax1.set_ylabel('残差')
            ax1.grid(True)
            
            # 2. 最大空泡率随时间的变化
            ax2.plot(self.history['time'], self.history['max_void_fraction'])
            ax2.set_title('最大空泡率随时间的变化')
            ax2.set_xlabel('时间 (s)')
            ax2.set_ylabel('最大空泡率')
            ax2.grid(True)
            
            # 3. 平均压力随时间的变化
            ax3.plot(self.history['time'], self.history['avg_pressure'])
            ax3.set_title('平均压力随时间的变化')
            ax3.set_xlabel('时间 (s)')
            ax3.set_ylabel('压力 (Pa)')
            ax3.grid(True)
            
            # 4. 最大温度随时间的变化
            ax4.plot(self.history['time'], self.history['max_temperature'])
            ax4.set_title('最大温度随时间的变化')
            ax4.set_xlabel('时间 (s)')
            ax4.set_ylabel('温度 (K)')
            ax4.grid(True)
            
            plt.tight_layout()
            plt.show()
            
            # 创建动画展示场分布随时间的变化
            self.animate_field_evolution()
            
        except Exception as e:
            print(f"绘制瞬态结果时出错: {str(e)}")

    def animate_field_evolution(self):
        """创建场分布随时间变化的动画并保存"""
        try:
            # 创建动画图形
            fig = plt.figure(figsize=(15, 5))
            
            # 创建三个子图：空泡率、温度和压力
            ax1 = fig.add_subplot(131, projection='3d')
            ax2 = fig.add_subplot(132, projection='3d')
            ax3 = fig.add_subplot(133, projection='3d')
            
            x = np.arange(self.n_axial)
            y = np.arange(self.n_channels)
            X, Y = np.meshgrid(x, y)
            
            # 确保有足够的历史数据
            n_frames = min(len(self.history['time']),
                      len(self.history['void_fraction_history']),
                      len(self.history['temperature_history']),
                      len(self.history['pressure_history']))
            
            if n_frames == 0:
                print("警告：没有可用的历史数据来创建动画")
                return
                
            def update(frame):
                # 清除当前帧
                ax1.clear()
                ax2.clear()
                ax3.clear()
                
                # 确保frame不会超出数组范围
                frame = min(frame, n_frames - 1)
                
                try:
                    # 绘制空泡率分布
                    surf1 = ax1.plot_surface(X, Y, 
                        self.history['void_fraction_history'][frame].T,
                        cmap=cm.viridis)
                    ax1.set_title(f'空泡率分布\n时间: {self.history["time"][frame]:.3f} s')
                    
                    # 绘制温度分布
                    surf2 = ax2.plot_surface(X, Y,
                        self.history['temperature_history'][frame].T,
                        cmap=cm.viridis)
                    ax2.set_title('温度分布 (K)')
                    
                    # 绘制压力分布
                    surf3 = ax3.plot_surface(X, Y,
                        self.history['pressure_history'][frame].T,
                        cmap=cm.viridis)
                    ax3.set_title('压力分布 (Pa)')
                    
                    # 设置视角
                    for ax in [ax1, ax2, ax3]:
                        ax.view_init(elev=20, azim=45)
                        ax.grid(True)
                        ax.set_xlabel('轴向位置')
                        ax.set_ylabel('子通道编号')
                        
                except IndexError as e:
                    print(f"警告：动画帧 {frame} 超出数据范围")
                    return
            
            # 创建动画
            anim = FuncAnimation(fig, update, 
                                frames=n_frames,
                                interval=200)
            
            # 保存动画
            try:
                anim.save(f'{self.results_folder}/evolution.mp4', 
                         writer='ffmpeg', fps=10)
                print(f"动画已保存到 {self.results_folder}/evolution.mp4")
            except Exception as e:
                print(f"保存动画失败: {str(e)}")
                print("尝试使用其他格式保存...")
                try:
                    anim.save(f'{self.results_folder}/evolution.gif', 
                             writer='pillow', fps=10)
                    print(f"动画已保存为GIF格式")
                except Exception as e:
                    print(f"保存GIF也失败: {str(e)}")
            
            plt.tight_layout()
            plt.show()
            
        except Exception as e:
            print(f"创建动画时出错: {str(e)}")

    def create_results_folder(self):
        """创建结果保存文件夹"""
        if not os.path.exists(self.results_folder):
            os.makedirs(self.results_folder)
            

    def solve_mass_momentum_equations(self):
        """求解质量和动量方程 - Numba优化版本"""
        # 使用Numba优化的网格点计算
        new_void, new_vl, new_vv = compute_grid_points_numba(
            self.n_axial, self.n_channels,
            self.void_fraction, self.prev_void_fraction,
            self.vapor_velocity, self.liquid_velocity,
            self.prev_liquid_velocity, self.prev_vapor_velocity,
            self.pressure, self.temperature,
            self.time_step, self.dz
        )
        
        # 更新结果
        self.void_fraction = new_void
        self.liquid_velocity = new_vl
        self.vapor_velocity = new_vv
        
    def solve_energy_equations(self):
        """求解能量方程，适应16子通道结构"""
        print("\n=== 开始求解能量方程 ===")
        try:
            dt = 1e-4  # 时间步长
            print(f"时间步长: {dt}")
            
            for i in range(1, self.n_axial):
                for j in range(self.n_channels):
                    print(f"\n处理子通道 {j+1} 在轴向位置 {i*self.dz:.3f} m")
                    
                    # 获取子通道类型和几何参数
                    channel_type = self.channel_types[j+1]
                    flow_area = self.channel_geometry[channel_type]['flow_area']
                    heated_perimeter = self.channel_geometry[channel_type]['heated_perimeter']
                    
                    print(f"子通道类型: {channel_type}")
                    print(f"流通面积: {flow_area:.6f} m²")
                    print(f"加热周长: {heated_perimeter:.6f} m")
                    
                    # 计算横向能量传递
                    net_energy_transfer = 0.0
                    for k in self.channel_connections[j+1]:
                        # 温度混合
                        w_prime = self.w_prime_ij[j, k-1]
                        dT = self.T_l[i, k-1] - self.T_l[i, j]
                        
                        print(f"与子通道 {k} 的温度差: {dT:.2f} K")
                        print(f"湍流混合系数: {w_prime:.4f}")
                        
                        # 获取流体物性
                        props_l = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
                        if props_l is not None:
                            net_energy_transfer += w_prime * props_l.get('specific_heat', 4200) * dT
                    
                    print(f"净横向能量传递: {net_energy_transfer:.2f} W")
                    
                    # 计算壁面换热
                    h_l = self.calculate_liquid_htc(i, j)  # 液相换热系数
                    h_v = self.calculate_vapor_htc(i, j)   # 气相换热系数
                    
                    print(f"液相换热系数: {h_l:.2f} W/m²·K")
                    print(f"气相换热系数: {h_v:.2f} W/m²·K")
                    
                    # 计算相变率
                    gamma = self.calculate_phase_change_rate(i, j)
                    print(f"相变率: {gamma:.4f} kg/m³·s")
                    
                    # 获取气相物性
                    props_v = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')
                    if props_v is None:
                        print("警告: 无法获取气相物性")
                        continue
                    
                    # 液相���量方程
                    dh_l = (
                        -self.liquid_velocity[i,j] * (self.h_l[i,j] - self.h_l[i-1,j]) / self.dz  # 对流
                        + h_l * heated_perimeter * (self.T_c[i,j] - self.T_l[i,j]) / (props_l['density'] * flow_area)  # 壁面换热
                        + net_energy_transfer / (props_l['density'] * flow_area)  # 横向混合
                        - gamma * self.h_l[i,j] / (props_l['density'] * (1 - self.void_fraction[i,j]))  # 相变
                    ) * dt
                    
                    print(f"液相焓变化: {dh_l:.2f} J/kg")
                    
                    # 气相能量方程
                    dh_v = (
                        -self.vapor_velocity[i,j] * (self.h_v[i,j] - self.h_v[i-1,j]) / self.dz  # 对流
                        + h_v * heated_perimeter * (self.T_c[i,j] - self.T_v[i,j]) / (props_v['density'] * flow_area)  # 壁面换热
                        + gamma * self.h_v[i,j] / (props_v['density'] * self.void_fraction[i,j])  # 相变
                    ) * dt
                    
                    print(f"气相焓变化: {dh_v:.2f} J/kg")
                    
                    # 更新焓值
                    self.h_l[i,j] += dh_l
                    self.h_v[i,j] += dh_v
                    
                    # 更新温度
                    self.T_l[i,j] = self.get_temperature_from_enthalpy(self.h_l[i,j], self.pressure[i,j], 'liquid')
                    self.T_v[i,j] = self.get_temperature_from_enthalpy(self.h_v[i,j], self.pressure[i,j], 'vapor')
                    
                    print(f"更新后的液相温度: {self.T_l[i,j]:.2f} K")
                    print(f"更新后的气相温度: {self.T_v[i,j]:.2f} K")
                    
        except Exception as e:
            print(f"求解能量方程时出错: {str(e)}")
            print(f"错误位置: {traceback.format_exc()}")
    
    def calculate_phase_change_rate(self, i, j):
        """
        计算相变率，考虑多种相变机制
        
        包括:
        1. 过冷沸腾
        2. 饱和沸腾
        3. 凝结
        """
        try:
            # 获取局部条件
            P = self.pressure[i,j]
            T_l = self.T_l[i,j]
            T_v = self.T_v[i,j]
            T_w = self.T_c[i,j]
            
            # 使用CoolProp获取精确的热力学参数
            T_sat = PropsSI('T', 'P', P, 'Q', 0, 'Water')
            h_fg = PropsSI('H', 'P', P, 'Q', 1, 'Water') - PropsSI('H', 'P', P, 'Q', 0, 'Water')
            
            # 获取流体物性
            props_l = self.get_fluid_properties(P, T_l, 'liquid')
            props_v = self.get_fluid_properties(P, T_v, 'vapor')
            
            if props_l is None or props_v is None:
                return 0.0
                
            # 1. 过冷沸腾
            if T_l < T_sat and T_w > T_sat:
                # 使用改进的过冷沸腾模型
                q_sc = self.calculate_subcooled_boiling_heat_flux(T_w, T_sat, T_l, P)
                gamma_sc = q_sc / h_fg
                
            # 2. 饱和沸腾
            elif T_l >= T_sat and T_w > T_sat:
                # 使用改进的核态沸腾模型
                q_nb = self.calculate_nucleate_boiling_heat_flux(T_w, T_sat, P)
                gamma_nb = q_nb / h_fg
                
            # 3. 凝结
            elif T_v > T_sat:
                # 使用改进的凝结模型
                q_cond = self.calculate_condensation_heat_flux(T_v, T_sat, P)
                gamma_cond = -q_cond / h_fg
                
            else:
                return 0.0
                
            # 考虑界面面积密度的影响
            a_i = 6 * self.void_fraction[i,j] / 0.001  # 假设气泡直径为1mm
            
            # 计算总的相变率
            if T_l < T_sat and T_w > T_sat:
                gamma = gamma_sc * a_i
            elif T_l >= T_sat and T_w > T_sat:
                gamma = gamma_nb * a_i
            else:
                gamma = gamma_cond * a_i
                
            return gamma
            
        except Exception as e:
            print(f"计算相变率时出错: {str(e)}")
            return 0.0

    def calculate_subcooled_boiling_heat_flux(self, T_w, T_sat, T_l, P):
        """改进的过冷沸腾换热模型"""
        try:
            # 获取流体物性
            props_l = self.get_fluid_properties(P, T_l, 'liquid')
            if props_l is None:
                return 0.0
                
            # 计算对流换热系数
            h_c = self.calculate_single_phase_htc(T_w, P, 'liquid')
            
            # 计算沸腾抑制因子
            S = max(0, min(1, (T_w - T_sat)/(T_sat - T_l)))
            
            # 计算核态沸腾换热系数
            h_nb = self.calculate_nucleate_boiling_htc(T_w, T_sat, P)
            
            # 组合换热系数
            q_sc = h_c * (T_w - T_l) + S * h_nb * (T_w - T_sat)
            
            return q_sc
            
        except Exception as e:
            print(f"计算过冷沸腾热流密度时出错: {str(e)}")
            return 0.0

    def calculate_nucleate_boiling_heat_flux(self, T_w, T_sat, P):
        """改进的核态沸腾换热模型"""
        try:
            # 使用改进的Rohsenow相关式
            props_l = self.get_fluid_properties(P, T_sat, 'liquid')
            props_v = self.get_fluid_properties(P, T_sat, 'vapor')
            
            if props_l is None or props_v is None:
                return 0.0
                
            h_fg = PropsSI('H', 'P', P, 'Q', 1, 'Water') - PropsSI('H', 'P', P, 'Q', 0, 'Water')
            sigma = 0.072  # 水的表面张力系数
            g = 9.81
            
            # 计算无量纲参数
            Ja = props_l['specific_heat'] * (T_w - T_sat) / h_fg
            Ra = (g * (props_l['density'] - props_v['density']) * 
                  props_l['density']**2 * h_fg * self.channel_diameter**3 /
                  (props_l['viscosity'] * props_l['thermal_conductivity'] * (T_w - T_sat)))
                  
            # 计算核态沸腾热流密度
            q_nb = props_l['thermal_conductivity'] * (T_w - T_sat) / self.channel_diameter * \
                   (Ja * Ra)**0.33
                   
            return q_nb
            
        except Exception as e:
            print(f"计算核态沸腾热流密度时出错: {str(e)}")
            return 0.0

    def calculate_condensation_heat_flux(self, T_v, T_sat, P):
        """改进的凝结换热模型"""
        try:
            # 使用Nusselt凝结理论
            props_l = self.get_fluid_properties(P, T_sat, 'liquid')
            props_v = self.get_fluid_properties(P, T_sat, 'vapor')
            
            if props_l is None or props_v is None:
                return 0.0
                
            h_fg = PropsSI('H', 'P', P, 'Q', 1, 'Water') - PropsSI('H', 'P', P, 'Q', 0, 'Water')
            g = 9.81
            
            # 计算凝结膜厚度
            delta = (4 * props_l['viscosity'] * props_l['thermal_conductivity'] * 
                    (T_v - T_sat) * self.channel_length /
                    (g * props_l['density'] * (props_l['density'] - props_v['density']) * h_fg))**0.25
                    
            # 计算凝结热流密度
            q_cond = props_l['thermal_conductivity'] * (T_v - T_sat) / delta
            
            return q_cond
            
        except Exception as e:
            print(f"计算凝结热流密度时出错: {str(e)}")
            return 0.0

    def calculate_turbulent_mixing(self, i, j, phase):
        """计算湍流混合和空泡漂移质量流量"""
        C_e = 0.01  # 涡流扩散率
        L_tur = self.channel_diameter  # 湍流混合长度
        
        mixing_mass = 0.0
        for k in self.channel_connections[j+1]:
            L_gap = self.get_gap_width(j+1, k)  # 获取间隙宽度
            
            if phase == 'liquid':
                rho_diff = self.get_fluid_properties(self.pressure[i,j], 
                    self.temperature[i,j], 'liquid')['density']
                alpha_diff = self.alpha_l[i,j] - self.alpha_l[i,k-1]
                alpha_eq_diff = self.calculate_equilibrium_void_fraction(i,j) - \
                               self.calculate_equilibrium_void_fraction(i,k-1)
            else:
                rho_diff = self.get_fluid_properties(self.pressure[i,j],
                    self.temperature[i,j], 'vapor')['density']
                alpha_diff = self.alpha_v[i,j] - self.alpha_v[i,k-1]
                alpha_eq_diff = (1-self.calculate_equilibrium_void_fraction(i,j)) - \
                               (1-self.calculate_equilibrium_void_fraction(i,k-1))
                
            mixing_mass += C_e * L_gap * L_tur * rho_diff * (alpha_diff - alpha_eq_diff)
            
        return mixing_mass

    def calculate_wall_friction(self, i, j, phase, direction='axial'):
        """计算壁面摩擦力"""
        print(f"\n=== 计算壁面摩擦力 (i={i}, j={j}, phase={phase}, direction={direction}) ===")
        if phase == 'liquid':
            props = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
            if direction == 'axial':
                velocity = self.liquid_velocity[i,j]
            else:
                velocity = self.v_l[i,j]
        else:
            props = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')
            if direction == 'axial':
                velocity = self.vapor_velocity[i,j]
            else:
                velocity = self.v_v[i,j]
                
        if props is None:
            print("警告: 无法获取流体物性")
            return 0.0
            
        print(f"流体速度: {velocity:.4f} m/s")
        print(f"流体密度: {props['density']:.2f} kg/m³")
        print(f"流体粘度: {props['viscosity']:.6f} Pa·s")
        
        friction = calculate_wall_friction_numba(
            velocity, self.channel_diameter,
            props['density'], props['viscosity']
        )
        
        print(f"计算得到的壁面摩擦力: {friction:.4f} N/m²")
        return friction

    def calculate_residual(self):
        """计算残差 - Numba优化版本"""
        return calculate_residual_numba(
            self.n_axial, self.n_channels,
            self.liquid_velocity, self.liquid_enthalpy,
            self.dz
        )

    def get_latent_heat(self, pressure):
        """获取给定压力下的汽化潜热"""
        try:
            h_f = PropsSI('H', 'P', pressure, 'Q', 0, 'Water')  # 饱和液焓
            h_g = PropsSI('H', 'P', pressure, 'Q', 1, 'Water')  # 饱和气焓
            return h_g - h_f
        except:
            return 2.257e6  # 返回1个大气压下的近似值

    def get_saturation_temperature(self, pressure):
        """
        计算给定压力下的饱和温度
        
        参数：
        pressure: float, 压力 (Pa)
        
        返回：
        float: 饱和温度 (K)
        """
        try:
            return PropsSI('T', 'P', pressure, 'Q', 0, 'Water')
        except Exception as e:
            print(f"警告：饱和温度计算失败 (P={pressure})")
            print(f"错误信息: {str(e)}")
            return 373.15  # 返回1个大气压下的近似值

    def calculate_single_phase_htc(self, i, j, phase):
        """
        计算单相对流换热系数
        
        参数：
        i, j: int, 网格位置
        phase: str, 'liquid' 或 'vapor'
        
        返回：
        float: 换热系数 (W/m²·K)
        """
        if phase == 'liquid':
            props = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
            velocity = self.liquid_velocity[i,j]
        else:
            props = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')
            velocity = self.vapor_velocity[i,j]
        
        if props is None:
            return 0.0
        
        Re = props['density'] * velocity * self.channel_diameter / props['viscosity']
        Pr = props['viscosity'] * props.get('specific_heat', 4200) / props.get('thermal_conductivity', 0.6)
        
        # Dittus-Boelter correlation
        if Re > 10000 and 0.7 <= Pr <= 160:
            Nu = 0.023 * Re**0.8 * Pr**0.4
        else:
            Nu = 4.36  # 层流或过渡区域的简化处理
        
        h = Nu * props.get('thermal_conductivity', 0.6) / self.channel_diameter
        return h

    def check_convergence(self, old_void, old_vl, old_vv, old_p, old_t):
        """检查收敛性"""
        max_error = max(
            np.max(np.abs(self.void_fraction - old_void)),
            np.max(np.abs(self.liquid_velocity - old_vl)),
            np.max(np.abs(self.vapor_velocity - old_vv)),
            np.max(np.abs(self.pressure - old_p)) / self.inlet_pressure,
            np.max(np.abs(self.temperature - old_t)) / self.inlet_temp
        )
        return max_error

    def setup_staggered_grid(self):
        """设置交错网格"""
        # 标量网格（压力、温度、空泡率等）
        self.scalar_grid = np.zeros((self.n_axial, self.n_channels))
        
        # 轴向速度网格（错开半个网格）
        self.axial_velocity_grid = np.zeros((self.n_axial + 1, self.n_channels))
        
        # 横向速度网格（错开半个网格）
        self.lateral_velocity_grid = np.zeros((self.n_axial, self.n_channels + 1))
        
        # 网格中心点坐标
        self.z_centers = np.linspace(0, self.channel_length, self.n_axial)
        self.z_faces = np.linspace(-self.dz/2, self.channel_length + self.dz/2, self.n_axial + 1)
        
    def build_preconditioner(self):
        """构建基于物理过程的预处理器"""
        n = self.n_axial * self.n_channels
        M = lil_matrix((n, n))
        
        # 构建预处理器矩阵（仅考虑轴向相邻单元）
        for i in range(self.n_axial):
            for j in range(self.n_channels):
                idx = i * self.n_channels + j
                
                # 中心单元
                M[idx, idx] = 1.0
                
                # 轴向相邻单元
                if i > 0:
                    M[idx, (i-1) * self.n_channels + j] = -0.5
                if i < self.n_axial - 1:
                    M[idx, (i+1) * self.n_channels + j] = -0.5
                    
        return M.tocsr()
        
    def solve_pkm(self):
        """使用Picard Krylov方法求解方程组"""
        converged = False
        picard_iter = 0
        
        while not converged and picard_iter < self.max_picard_iterations:
            # 1. 计算热通量
            self.calculate_heat_flux()
            
            # 2. 求解动量方程
            self.solve_momentum_equations()
            
            # 3. 构建质量-能量方程组
            A, b = self.build_mass_energy_system()
            
            # 4. 构建预处理器
            M = self.build_preconditioner()
            
            # 5. 使用PGMRES求解线性方程组
            try:
                x, info = gmres(A, b, M=M, tol=self.krylov_tolerance, 
                              maxiter=self.max_krylov_iterations)
                
                if info > 0:
                    print(f"警告：GMRES在{info}次迭代后未收敛")
                    
                # 更新解
                self.update_solution(x)
                
                # 计算收敛性
                residual = self.calculate_residual()
                converged = residual < self.picard_tolerance
                
            except Exception as e:
                print(f"错误：GMRES求解失败 - {str(e)}")
                # 如果求解失败，减小时间步长重试
                self.time_step *= 0.5
                if self.time_step < 1e-6:
                    raise Exception("时间步长过小，计算终止")
                continue
                
            picard_iter += 1
            
        return converged, picard_iter
        
    def build_mass_energy_system(self):
        """构建质量-能量方程组的系数矩阵和右端向量"""
        n = self.n_axial * self.n_channels
        A = lil_matrix((n, n))
        b = np.zeros(n)
        
        for i in range(self.n_axial):
            for j in range(self.n_channels):
                idx = i * self.n_channels + j
                
                # 1. 质量守恒方程
                # 轴向质量流量
                if i > 0:
                    A[idx, idx-self.n_channels] = -self.liquid_velocity[i-1,j]
                A[idx, idx] = self.liquid_velocity[i,j]
                if i < self.n_axial - 1:
                    A[idx, idx+self.n_channels] = -self.liquid_velocity[i+1,j]
                
                # 横向质量流量
                for k in self.channel_connections[j+1]:
                    k_idx = i * self.n_channels + (k-1)
                    w_jk = self.calculate_cross_flow(i, j+1, k)
                    A[idx, k_idx] = w_jk
                    
                # 2. 能量守恒方程
                # 轴向能量输运
                if i > 0:
                    A[idx, idx-self.n_channels] += -self.liquid_velocity[i-1,j] * self.liquid_enthalpy[i-1,j]
                A[idx, idx] += self.liquid_velocity[i,j] * self.liquid_enthalpy[i,j]
                if i < self.n_axial - 1:
                    A[idx, idx+self.n_channels] += -self.liquid_velocity[i+1,j] * self.liquid_enthalpy[i+1,j]
                
                # 加热项
                q_wall = self.calculate_wall_heat_flux(i, j, 'liquid')
                b[idx] = q_wall
                
        return A.tocsr(), b
        
    def solve_momentum_equations(self):
        """求解动量方程"""
        # 轴向动量方程
        for i in range(1, self.n_axial):
            for j in range(self.n_channels):
                # 获取交错网格上的速度
                v_up = self.axial_velocity_grid[i-1, j]
                v_center = self.axial_velocity_grid[i, j]
                v_down = self.axial_velocity_grid[i+1, j] if i < self.n_axial-1 else v_center
                
                # 计算压力梯度
                dp_dz = (self.pressure[i,j] - self.pressure[i-1,j]) / self.dz
                
                # 计算摩擦力
                f_wall = self.calculate_wall_friction(i, j, 'liquid')
                
                # 更新轴向速度
                self.axial_velocity_grid[i,j] = v_center - self.time_step * (
                    v_center * (v_down - v_up)/(2*self.dz) +  # 对流项
                    dp_dz/self.get_fluid_properties(self.pressure[i,j], 
                        self.temperature[i,j], 'liquid')['density'] +  # 压力项
                    f_wall +  # 摩擦力
                    9.81  # 重力项
                )
        
        # 横向动量方程
        for i in range(self.n_axial):
            for j in range(self.n_channels-1):
                # 计算横向压力梯度
                dp_dy = self.calculate_lateral_pressure_gradient(i, j)
                
                # 计算横向摩擦力
                f_lat = self.calculate_wall_friction(i, j, 'liquid', 'lateral')
                
                # 更新横向速度
                self.lateral_velocity_grid[i,j] = self.lateral_velocity_grid[i,j] - self.time_step * (
                    dp_dy/self.get_fluid_properties(self.pressure[i,j], 
                        self.temperature[i,j], 'liquid')['density'] +
                    f_lat
                )
                
    def update_solution(self, x):
        """更新解向量"""
        for i in range(self.n_axial):
            for j in range(self.n_channels):
                idx = i * self.n_channels + j
                
                # 更新压力
                self.pressure[i,j] = x[idx]
                
                # 更新温度（通过压力和焓计算）
                h = x[idx + self.n_axial * self.n_channels]
                self.temperature[i,j] = self.get_temperature_from_enthalpy(
                    h, self.pressure[i,j], 'liquid')
                
    def calculate_heat_flux(self):
        """计算热通量"""
        # 计算热通量
        # ...

    def calculate_friction_factor(self, velocity, viscosity, flow_area, wetted_perimeter):
        """计算摩擦因子"""
        # 计算水力直径
        hydraulic_diameter = 4 * flow_area / wetted_perimeter
        
        # 计算雷诺数
        reynolds = abs(velocity) * hydraulic_diameter / viscosity
        
        # 层流区
        if reynolds < 2300:
            return 64 / reynolds
        # 湍流区
        else:
            return 0.316 * reynolds**(-0.25)

    def debug_output(self, time_step):
        """输出调试信息，用于验证16子通道计算"""
        try:
            print(f"\n====== 时间步 {time_step} 调试信息 ======")
            
            # 1. 检查质量守恒
            inlet_mass_flow = 0
            outlet_mass_flow = 0
            for j in range(self.n_channels):
                channel_type = self.channel_types[j+1]
                flow_area = self.channel_geometry[channel_type]['flow_area']
                
                # 入口质量流量
                rho_in = self.get_fluid_properties(self.pressure[0,j], self.T_l[0,j], 'liquid')['density']
                inlet_mass_flow += rho_in * self.liquid_velocity[0,j] * flow_area
                
                # 出口质量流量
                rho_l_out = self.get_fluid_properties(self.pressure[-1,j], self.T_l[-1,j], 'liquid')['density']
                rho_v_out = self.get_fluid_properties(self.pressure[-1,j], self.T_v[-1,j], 'vapor')['density']
                outlet_mass_flow += (rho_l_out * (1-self.void_fraction[-1,j]) * self.liquid_velocity[-1,j] + 
                                    rho_v_out * self.void_fraction[-1,j] * self.vapor_velocity[-1,j]) * flow_area
            
            mass_error = abs(inlet_mass_flow - outlet_mass_flow) / inlet_mass_flow
            print(f"质量守恒误差: {mass_error:.6f}")
            
            # 2. 检查能量守恒
            total_power_in = 0
            total_power_out = 0
            total_heat_added = 0
            
            for j in range(self.n_channels):
                channel_type = self.channel_types[j+1]
                flow_area = self.channel_geometry[channel_type]['flow_area']
                heated_perimeter = self.channel_geometry[channel_type]['heated_perimeter']
                
                # 入口能量
                h_in = self.h_l[0,j]
                total_power_in += inlet_mass_flow * h_in / self.n_channels
                
                # 出口能量
                h_l_out = self.h_l[-1,j]
                h_v_out = self.h_v[-1,j]
                total_power_out += outlet_mass_flow * (
                    (1-self.void_fraction[-1,j]) * h_l_out + 
                    self.void_fraction[-1,j] * h_v_out
                ) / self.n_channels
                
                # 加热功率
                for i in range(self.n_axial):
                    total_heat_added += self.q_wl[i,j] * self.dz
            
            energy_error = abs(total_power_out - total_power_in - total_heat_added) / total_heat_added
            print(f"能量守恒误差: {energy_error:.6f}")
            
            # 3. 检查各子通道类型的平均值
            for type_name in ['corner', 'edge', 'center']:
                void_sum = 0
                temp_sum = 0
                count = 0
                
                for j in range(self.n_channels):
                    if self.channel_types[j+1] == type_name:
                        void_sum += np.mean(self.void_fraction[:,j])
                        temp_sum += np.mean(self.T_l[:,j])
                        count += 1
                
                if count > 0:
                    print(f"\n{type_name}型子通道平均值:")
                    print(f"平均空泡率: {void_sum/count:.4f}")
                    print(f"平均温度: {temp_sum/count:.2f} K")
            
            # 4. 检查数值稳定性
            max_courant = 0
            for j in range(self.n_channels):
                channel_type = self.channel_types[j+1]
                flow_area = self.channel_geometry[channel_type]['flow_area']
                
                # 计算Courant数
                max_velocity = max(np.max(abs(self.liquid_velocity[:,j])), 
                                 np.max(abs(self.vapor_velocity[:,j])))
                courant = max_velocity * self.dt / self.dz
                max_courant = max(max_courant, courant)
            
            print(f"\n最��Courant数: {max_courant:.4f}")
            if max_courant > 1:
                print("警告：Courant数大于1，可能存在数值不稳定")
                
        except Exception as e:
            print(f"尝试输出时出错: {str(e)}")

    def verify_channel_connections(self):
        """验证子通道连接的正确性"""
        try:
            print("\n====== 验证子通道连接 ======")
            
            # 检查连接的对称性
            for j in range(1, self.n_channels + 1):
                for k in self.channel_connections[j]:
                    if j not in self.channel_connections[k]:
                        print(f"错误：子通道{j}和{k}的连接不对称")
                        
            # 检查连接的物理合理性
            for j in range(1, self.n_channels + 1):
                connections = self.channel_connections[j]
                channel_type = self.channel_types[j]
                
                # 检查连接数量
                if channel_type == 'corner' and len(connections) != 2:
                    print(f"警告：角通道{j}的连接数量不正确")
                elif channel_type == 'edge' and len(connections) != 3:
                    print(f"警告：边通道{j}的连接数量不正确")
                elif channel_type == 'center' and len(connections) != 4:
                    print(f"警告：中心通道{j}的连接数量不正确")
                    
            print("子通道连接验证完成")
            
        except Exception as e:
            print(f"验证子通道连接时出错: {str(e)}")

    def check_numerical_stability(self):
        """检查数值计算的稳定性"""
        try:
            print("\n====== 检查数值稳定性 ======")
            
            # 检查物理量的范围
            if np.any(self.void_fraction < 0) or np.any(self.void_fraction > 1):
                print("警告：空泡率超出物理范围")
                
            if np.any(self.pressure < 0):
                print("警告：出现负压力")
                
            if np.any(self.T_l < 273) or np.any(self.T_v < 273):
                print("警告：温度低于273K")
                
            # 检查数值梯度
            max_void_gradient = np.max(np.abs(np.diff(self.void_fraction, axis=0)))
            max_pressure_gradient = np.max(np.abs(np.diff(self.pressure, axis=0)))
            
            print(f"最大空泡率梯度: {max_void_gradient:.6f}")
            print(f"最大压力梯度: {max_pressure_gradient:.2f} Pa/m")
            
            if max_void_gradient > 0.1:
                print("警告：空泡率梯度过大")
            if max_pressure_gradient > 1e5:
                print("警告：压力梯度过大")
                
        except Exception as e:
            print(f"检查数值稳定性时出错: {str(e)}")

    def calculate_liquid_htc(self, i, j):
        """
        计算液相换热系数
        
        使用Dittus-Boelter关联式计算单相强制对流换热系数。
        考虑不同类型子通道的几何特征。
        
        参数：
            i (int): 轴向位置索引
            j (int): 子通道编号
            
        返回：
            float: 液相换热系数 (W/m²·K)
            
        注意：
            - 使用Dittus-Boelter关联式: Nu = 0.023 * Re^0.8 * Pr^0.4
            - 考虑不同类型子通道的水力直径
            - 包含异常处理机制
        """
        try:
            # 获取子通道类型和几何参数
            channel_type = self.channel_types[j]
            flow_area = self.channel_geometry[channel_type]['flow_area']
            wetted_perimeter = self.channel_geometry[channel_type]['wetted_perimeter']
            
            # 计算水力直径
            d_h = 4.0 * flow_area / wetted_perimeter
            
            # 获取流体物性
            p = self.pressure[i, j-1]
            t = self.temperature[i, j-1]
            props = self.get_fluid_properties(p, t, 'liquid')
            
            if props is None:
                return 0.0
                
            rho = props['density']
            mu = props['viscosity']
            
            # 计算雷诺数
            v = self.liquid_velocity[i, j-1]
            re = rho * abs(v) * d_h / mu
            
            # 计算普朗特数 (假设为常数)
            pr = 1.0
            
            # 使用Dittus-Boelter关联式
            if re > 2300:  # 湍流
                nu = 0.023 * re**0.8 * pr**0.4
            else:  # 层流
                nu = 3.66
                
            # 计算导热系数 (假设为常数)
            k = 0.6  # W/m·K
            
            # 计算换热系数
            h = nu * k / d_h
            
            return h
            
        except Exception as e:
            print(f"计算液相换热系数时出错: {str(e)}")
            return 0.0
            
    def calculate_vapor_htc(self, i, j):
        """
        计算气相换热系数
        
        使用改进的Dittus-Boelter关联式计算气相的换热系数。
        考虑气相的特殊流动特性。
        
        参数：
            i (int): 轴向位置索引
            j (int): 子通道编号
            
        返回：
            float: 气相换热系数 (W/m²·K)
            
        注意：
            - 使用修正的关联式: Nu = 0.018 * Re^0.8 * Pr^0.4
            - 系数0.018考虑了气相流动的特殊性
            - 包含对低雷诺数情况的处理
        """
        try:
            # 获取子通道几何参数
            channel_type = self.channel_types[j]
            flow_area = self.channel_geometry[channel_type]['flow_area']
            wetted_perimeter = self.channel_geometry[channel_type]['wetted_perimeter']
            
            # 计算水力直径
            d_h = 4.0 * flow_area / wetted_perimeter
            
            # 获取流体物性
            p = self.pressure[i, j-1]
            t = self.temperature[i, j-1]
            props = self.get_fluid_properties(p, t, 'vapor')
            
            if props is None:
                return 0.0
                
            rho = props['density']
            mu = props['viscosity']
            
            # 计算雷诺数
            v = self.vapor_velocity[i, j-1]
            re = rho * abs(v) * d_h / mu
            
            # 计算普朗特数 (假设为常数)
            pr = 1.0
            
            # 使用修正的关联式
            if re > 2300:  # 湍流
                nu = 0.018 * re**0.8 * pr**0.4  # 系数改为0.018
            else:  # 层流
                nu = 3.66
                
            # 计算导热系数 (假设为常数)
            k = 0.025  # W/m·K (气相导热系数较小)
            
            # 计算换热系数
            h = nu * k / d_h
            
            return h
            
        except Exception as e:
            print(f"计算气相换热系数时出错: {str(e)}")
            return 0.0
            
    def calculate_wall_temperature(self, i, j):
        """
        计算壁面温度
        
        基于热流密度和换热系数计算壁面温度。
        考虑两相流动的换热特性。
        
        参数：
            i (int): 轴向位置索引
            j (int): 子通道编号
            
        返回：
            float: 壁面温度 (K)
            
        注意：
            - 考虑液相和气相的共同换热效果
            - 使用加权平均的方法处理两相换热
            - 包含对极限情况的处理
        """
        try:
            # 获取流体温度
            t_f = self.temperature[i, j-1]
            
            # 获取热流密度
            q_wl = self.q_wl[i, j-1]  # 液相热流密度
            q_wv = self.q_wv[i, j-1]  # 气相热流密度
            
            # 计算换热系数
            h_l = self.calculate_liquid_htc(i, j)
            h_v = self.calculate_vapor_htc(i, j)
            
            # 获取空泡率
            alpha = self.void_fraction[i, j-1]
            
            # 计算总的换热系数
            h_total = (1 - alpha) * h_l + alpha * h_v
            
            # 计算总的热流密度
            q_total = q_wl + q_wv
            
            # 计算壁面温度
            if h_total > 0:
                t_w = t_f + q_total / h_total
            else:
                t_w = t_f + 10  # 默认温差10K
                
            return t_w
            
        except Exception as e:
            print(f"计算壁面温度时出错: {str(e)}")
            return self.temperature[i, j-1] + 10  # 返回一个安全的估计值

    def calculate_condensation_heat_flux(self, T_w, T_sat, P):
        """计算冷凝热流密度"""
        try:
            # 获取气相物性
            props_v = self.get_fluid_properties(P, T_sat, 'vapor')
            if props_v is None or 'thermal_conductivity' not in props_v:
                print(f"警告：流体物性数据不完整 (P={P}, T_sat={T_sat})")
                return 0.0
            
            # 计算冷凝换热系数
            h_cond = self.calculate_single_phase_htc(T_w, P, 'vapor')  # 假设使用气相的单相换热系数
            
            # 计算冷凝热流密度
            q_cond = h_cond * (T_w - T_sat)
            
            return max(0, q_cond)  # 确保热流密度非负
            
        except Exception as e:
            print(f"计算冷凝热流密度时出错: {str(e)}")
            return 0.0

    def calculate_boiling_heat_flux(self, T_w, T_sat, P):
        """计算汽化热流密度"""
        try:
            # 获取液相物性
            props_l = self.get_fluid_properties(P, T_sat, 'liquid')
            if props_l is None or 'thermal_conductivity' not in props_l:
                print(f"警告：流体物性数据不完整 (P={P}, T_sat={T_sat})")
                return 0.0
            
            # 计算汽化换热系数
            h_boil = self.calculate_single_phase_htc(T_w, P, 'liquid')  # 假设使用液相的单相换热系数
            
            # 计算汽化热流密度
            q_boiling = h_boil * (T_w - T_sat)
            
            return max(0, q_boiling)  # 确保热流密度非负
            
        except Exception as e:
            print(f"计算汽化热流密度时出错: {str(e)}")
            return 0.0

    def get_temperature_from_enthalpy(self, enthalpy, pressure, phase):
        """
        根据给定的焓和压力计算温度
        
        参数：
            enthalpy (float): 比焓 (J/kg)
            pressure (float): 压力 (Pa)
            phase (str): 相态 ('liquid' 或 'vapor')
        
        返回：
            float: 温度 (K)
        """
        try:
            from CoolProp.CoolProp import PropsSI

            # 检查焓和压力的有效性
            if enthalpy <= 0 or pressure <= 0:
                raise ValueError("焓和压力必须为正值")

            # 使用PropsSI检查给定压力下的饱和焓
            T_sat = self.get_saturation_temperature(pressure)
            h_sat = PropsSI('H', 'P', pressure, 'Q', 0, 'Water')  # 饱和液体焓

            if enthalpy < h_sat:
                print(f"警告：给定焓 {enthalpy} 小于饱和液体焓 {h_sat}，可能不合理。")
                return T_sat  # 返回饱和温度

            def objective(T):
                if phase == 'liquid':
                    return PropsSI('H', 'P', pressure, 'T', T, 'Water') - enthalpy
                elif phase == 'vapor':
                    return PropsSI('H', 'P', pressure, 'T', T, 'Water') - enthalpy
                else:
                    raise ValueError("无效的相态参数，需为 'liquid' 或 'vapor'")

            # 初始猜测温度
            T_initial = self.get_saturation_temperature(pressure)
            T_solution = opt.newton(objective, T_initial)
            return T_solution

        except Exception as e:
            print(f"警告：根据焓和压力计算温度失败 (H={enthalpy}, P={pressure})")
            print(f"错误信息: {str(e)}")
            return 373.15  # 返回默认温度

    def calculate_interface_tension(self, i, j):
        """计算界面张力效应"""
        try:
            # 获取局部条件
            alpha = self.void_fraction[i,j]
            d_bubble = 0.001  # 假设的气泡直径 (m)
            
            # 计算界面曲率
            if i > 0 and i < self.n_axial-1:
                d2alpha_dz2 = (self.void_fraction[i+1,j] - 2*self.void_fraction[i,j] + 
                             self.void_fraction[i-1,j]) / (self.dz**2)
            else:
                d2alpha_dz2 = 0
                
            # 计算界面张力
            F_st = -self.surface_tension * d2alpha_dz2
            
            # 考虑壁面接触角的影响
            theta = np.radians(self.contact_angle)
            F_st *= np.cos(theta)
            
            return F_st
            
        except Exception as e:
            print(f"计算界面张力时出错: {str(e)}")
            return 0.0

    def solve_k_epsilon(self, i, j):
        """求解k-ε湍流模型方程"""
        try:
            # 获取局部流动参数
            u = self.liquid_velocity[i,j]
            rho = self.get_fluid_properties(self.pressure[i,j], 
                                          self.temperature[i,j], 'liquid')['density']
            mu = self.get_fluid_properties(self.pressure[i,j], 
                                         self.temperature[i,j], 'liquid')['viscosity']
            
            # 计算速度梯度
            if i > 0 and i < self.n_axial-1:
                du_dy = (self.liquid_velocity[i,j+1] - self.liquid_velocity[i,j-1]) / (2*self.dz)
                
                # 计算湍流产生项
                P_k = mu * (du_dy**2)
                
                # 求解k方程
                self.k[i,j] += self.dt * (
                    P_k - self.epsilon[i,j] +
                    self.calculate_k_diffusion(i, j)
                )
                
                # 求解ε方程
                self.epsilon[i,j] += self.dt * (
                    self.C_1 * self.epsilon[i,j] * P_k / self.k[i,j] -
                    self.C_2 * self.epsilon[i,j]**2 / self.k[i,j] +
                    self.calculate_epsilon_diffusion(i, j)
                )
                
                # 计算湍流粘度
                mu_t = self.C_mu * rho * self.k[i,j]**2 / self.epsilon[i,j]
                
                return mu_t
                
        except Exception as e:
            print(f"求解k-ε模型时出错: {str(e)}")
            return 0.0

    def solve_heat_coupling(self, i, j):
        """求解热耦合效应"""
        try:
            # 计算燃料到包壳的导热
            q_fuel_clad = self.k_fuel * (self.T_fuel[i,j] - self.T_clad[i,j]) / self.gap_conductance
            
            # 计算包壳到冷却剂的导热
            q_clad_cool = self.k_clad * (self.T_clad[i,j] - self.temperature[i,j])
            
            # 更新温度场
            self.T_fuel[i,j] += self.dt * (
                self.calculate_fuel_heat_generation(i, j) - q_fuel_clad
            ) / (self.rho_fuel * self.cp_fuel)
            
            self.T_clad[i,j] += self.dt * (
                q_fuel_clad - q_clad_cool
            ) / (self.rho_clad * self.cp_clad)
            
            return q_clad_cool
            
        except Exception as e:  # 添加异常处理
            print(f"求解热耦合效应时出错: {str(e)}")
            return 0.0  # 返���一个默认值以防出错

    def update_values(self, prev_void_fraction, prev_liquid_velocity, prev_vapor_velocity, time_step):
        """更新值"""
        new_void = np.zeros_like(prev_void_fraction)
        new_vl = np.zeros_like(prev_liquid_velocity)
        new_vv = np.zeros_like(prev_vapor_velocity)

        for i in range(new_void.shape[0]):
            for j in range(new_void.shape[1]):
                d_alpha_dx = 0
                d_vl_dx = 0
                d_vv_dx = 0
                dp_dx = 0

                # 更新值
                new_void[i,j] = prev_void_fraction[i,j] - time_step * (
                    vapor_velocity[i,j] * d_alpha_dx
                )

                new_vl[i,j] = prev_liquid_velocity[i,j] - time_step * (
                    liquid_velocity[i,j] * d_vl_dx +
                    dp_dx +
                    9.81
                )

                new_vv[i,j] = prev_vapor_velocity[i,j] - time_step * (
                    vapor_velocity[i,j] * d_vv_dx +
                    dp_dx +
                    9.81
                )

                # 限制物理范围
                new_void[i,j] = min(max(new_void[i,j], 0.0), 0.9999)
                new_vl[i,j] = max(new_vl[i,j], 0.0)
                new_vv[i,j] = max(new_vv[i,j], 0.0)

        return new_void, new_vl, new_vv
            
        

@njit(parallel=True)
def compute_energy_points_numba(n_axial, n_channels, temperature, prev_temperature,
                              liquid_velocity, pressure, time_step, dz):
    """使用Numba优化的能量点计算"""
    new_temp = np.zeros_like(temperature)
    
    for i in prange(1, n_axial):
        for j in prange(n_channels):
            # 计算空间导数项
            if i > 0:
                d_T_dx = (temperature[i,j] - temperature[i-1,j]) / dz
            else:
                d_T_dx = 0
            
            # 更新温度
            new_temp[i,j] = prev_temperature[i,j] - time_step * (
                liquid_velocity[i,j] * d_T_dx
            )
            
    return new_temp

@njit
def calculate_wall_friction_numba(velocity, diameter, density, viscosity):
    """使用Numba优化的壁面摩擦力计算"""
    # 计算Reynolds数
    Re = density * abs(velocity) * diameter / viscosity
    
    # 计算摩擦系数
    if Re < 2300:  # 层流
        f = 64 / Re
    else:  # 湍流
        f = 0.316 * Re**(-0.25)
        
    # 计算壁面摩擦力
    tau_w = 0.5 * f * density * velocity * abs(velocity)
    
    return tau_w

@njit(parallel=True)
def calculate_residual_numba(n_axial, n_channels, liquid_velocity, liquid_enthalpy, dz):
    """使用Numba优化的残差计算"""
    mass_residual = np.zeros((n_axial, n_channels))
    energy_residual = np.zeros((n_axial, n_channels))
    
    for i in prange(1, n_axial):
        for j in prange(n_channels):
            # 轴向质量流量差
            dm_dz = (liquid_velocity[i,j] - liquid_velocity[i-1,j]) / dz
            mass_residual[i,j] = dm_dz
            
            # 轴向能量输运
            dh_dz = (liquid_enthalpy[i,j] - liquid_enthalpy[i-1,j]) / dz
            energy_residual[i,j] = dh_dz
            
    return max(np.max(np.abs(mass_residual)), np.max(np.abs(energy_residual)))

@njit(parallel=True)
def solve_two_fluid_equations_numba(n_axial, n_channels, 
                                  alpha_l, alpha_v, rho_l, rho_v,
                                  u_l, u_v, v_l, v_v, w_l, w_v,
                                  P, T_l, T_v, T_c, h_l, h_v,
                                  dx, dy, dz, dt):
    """使用Numba优化的双流体模型求解器"""
    # 计算CFL时间步长
    dt_cfl = calculate_cfl_timestep(u_l, u_v, v_l, v_v, dx, dy, dz)
    
    # 使用较小的时间步长
    dt = min(dt, dt_cfl)
    
    # 初始化新的场值
    new_alpha_l = np.zeros_like(alpha_l)
    new_alpha_v = np.zeros_like(alpha_v)
    new_u_l = np.zeros_like(u_l)
    new_u_v = np.zeros_like(u_v)
    new_v_l = np.zeros_like(v_l)
    new_v_v = np.zeros_like(v_v)
    new_h_l = np.zeros_like(h_l)
    new_h_v = np.zeros_like(h_v)
    
    for i in prange(1, n_axial-1):
        for j in prange(1, n_channels-1):
            # ====== 1. 计算相变质量流量 ======
            T_sat = calculate_saturation_temperature_numba(P[i,j])
            
            # 过冷沸腾
            if T_l[i,j] < T_sat and T_c[i,j] > T_sat:
                h_sc = calculate_subcooled_boiling_htc_numba(T_c[i,j], T_sat, T_l[i,j], P[i,j])
                Q_sc = h_sc * (T_c[i,j] - T_l[i,j])
                gamma_sc = Q_sc / calculate_latent_heat_numba(P[i,j])
            else:
                gamma_sc = 0.0
                
            # 饱和核态沸腾
            if T_l[i,j] >= T_sat and T_c[i,j] > T_sat:
                h_sh = calculate_nucleate_boiling_htc_numba(T_c[i,j], T_sat, P[i,j])
                Q_sh = h_sh * (T_c[i,j] - T_sat)
                gamma_sh = Q_sh / calculate_latent_heat_numba(P[i,j])
            else:
                gamma_sh = 0.0
                
            # 净相变质量流量
            gamma = gamma_sh - gamma_sc
            
            # ====== 2. 计算湍流混合项 ======
            # 液相湍流混合
            T_ml = calculate_turbulent_mixing_numba(alpha_l[i,j], alpha_l[i,j+1], 
                                                  rho_l[i,j], dy)
            # 气相湍流混合
            T_mv = calculate_turbulent_mixing_numba(alpha_v[i,j], alpha_v[i,j+1], 
                                                  rho_v[i,j], dy)
                    
            # ====== 3. 计算界面力 ======
            F_lv = calculate_interfacial_force_numba(alpha_l[i,j], rho_l[i,j],
                                                   u_l[i,j], u_v[i,j],
                                                   v_l[i,j], v_v[i,j])
            
            # ====== 4. 液相质量守恒方程 ======
            d_alpha_l_dt = -(1/rho_l[i,j]) * (
                # 轴向质量通量导数
                (alpha_l[i+1,j]*rho_l[i+1,j]*u_l[i+1,j] - 
                 alpha_l[i-1,j]*rho_l[i-1,j]*u_l[i-1,j])/(2*dz) +
                # 横向质量通量导数
                (alpha_l[i,j+1]*rho_l[i,j+1]*v_l[i,j+1] - 
                 alpha_l[i,j-1]*rho_l[i,j-1]*v_l[i,j-1])/(2*dy)
            ) - gamma/rho_l[i,j] + T_ml
            
            # ====== 5. 气相质量守恒方程 ======
            d_alpha_v_dt = -(1/rho_v[i,j]) * (
                # 轴向���量通量导数
                (alpha_v[i+1,j]*rho_v[i+1,j]*u_v[i+1,j] - 
                 alpha_v[i-1,j]*rho_v[i-1,j]*u_v[i-1,j])/(2*dz) +
                # 横向质量通量导数
                (alpha_v[i,j+1]*rho_v[i,j+1]*v_v[i,j+1] - 
                 alpha_v[i,j-1]*rho_v[i,j-1]*v_v[i,j-1])/(2*dy)
            ) + gamma/rho_v[i,j] + T_mv
            
            # ====== 6. 液相轴向动量方程 ======
            # 壁面摩擦力
            F_wl = calculate_wall_friction_numba(u_l[i,j], 0.01, rho_l[i,j], 1e-3)
            
            d_ul_dt = -(
                # 对流项
                u_l[i,j] * (u_l[i+1,j] - u_l[i-1,j])/(2*dz) +
                v_l[i,j] * (u_l[i,j+1] - u_l[i,j-1])/(2*dy) +
                # 压力梯度
                (1/rho_l[i,j]) * (P[i+1,j] - P[i-1,j])/(2*dz) +
                # 重力项
                9.81 +
                # 壁面摩擦力
                F_wl/rho_l[i,j] -
                # 界面力
                F_lv/rho_l[i,j] -
                # 相变动量
                gamma * u_l[i,j]/rho_l[i,j]
            )
            
            # ====== 7. 气相轴向动量方程 ======
            # 壁面摩擦力
            F_wv = calculate_wall_friction_numba(u_v[i,j], 0.01, rho_v[i,j], 1e-5)
            
            d_uv_dt = -(
                # 对流项
                u_v[i,j] * (u_v[i+1,j] - u_v[i-1,j])/(2*dz) +
                v_v[i,j] * (u_v[i,j+1] - u_v[i,j-1])/(2*dy) +
                # 压力梯度
                (1/rho_v[i,j]) * (P[i+1,j] - P[i-1,j])/(2*dz) +
                # 重力项
                9.81 +
                # 壁面摩擦力
                F_wv/rho_v[i,j] +
                # 界面力
                F_lv/rho_v[i,j] +
                # 相变动量
                gamma * u_v[i,j]/rho_v[i,j]
            )
            
            # ====== 8. 液相横向动量方程 ======
            d_vl_dt = -(
                # 对流项
                u_l[i,j] * (v_l[i+1,j] - v_l[i-1,j])/(2*dz) +
                v_l[i,j] * (v_l[i,j+1] - v_l[i,j-1])/(2*dy) -
                # 界面力
                F_lv/rho_l[i,j] -
                # 相变动量
                gamma * v_l[i,j]/rho_l[i,j]
            )
            
            # ====== 9. 气相横向动量方程 ======
            d_vv_dt = -(
                # 对流项
                u_v[i,j] * (v_v[i+1,j] - v_v[i-1,j])/(2*dz) +
                v_v[i,j] * (v_v[i,j+1] - v_v[i,j-1])/(2*dy) +
                # 压力梯度
                (1/rho_v[i,j]) * (P[i,j+1] - P[i,j-1])/(2*dy) +
                # 界面力
                F_lv/rho_v[i,j] +
                # 相变动量
                gamma * v_v[i,j]/rho_v[i,j]
            )
            
            # ====== 10. 液相能量方程 - 添加燃料热生成和径向传导 ======
            # 壁面换热
            q_wl = calculate_wall_heat_flux_numba(T_c[i,j], T_l[i,j], h_l[i,j])
            
            # 计算燃料热生成
            q_fuel = calculate_fuel_heat_generation_numba(i, n_axial, dz)
            
            # 计算径向热传导
            q_rad = calculate_radial_heat_conduction_numba(T_c[i,j], T_l[i,j], dx)
            
            d_hl_dt = -(
                # 对流项
                u_l[i,j] * (h_l[i+1,j] - h_l[i-1,j])/(2*dz) +
                v_l[i,j] * (h_l[i,j+1] - h_l[i,j-1])/(2*dy) +
                # 压力功
                alpha_l[i,j] * (P[i,j] - P[i-1,j])/dt +
                # 壁面换热
                q_wl/(rho_l[i,j] * alpha_l[i,j]) +
                # 相变热量
                gamma * h_l[i,j]/(rho_l[i,j] * alpha_l[i,j]) +
                # 添加燃料热生成和径向传导
                (q_fuel + q_rad)/(rho_l[i,j] * alpha_l[i,j])
            )
            
            # ====== 11. 气相能量方程 - 添加燃料热生成和径向传导 ======
            # 壁面换热
            q_wv = calculate_wall_heat_flux_numba(T_c[i,j], T_v[i,j], h_v[i,j])
            
            d_hv_dt = -(
                # 对流项
                u_v[i,j] * (h_v[i+1,j] - h_v[i-1,j])/(2*dz) +
                v_v[i,j] * (h_v[i,j+1] - h_v[i,j-1])/(2*dy) +
                # 压力功
                alpha_v[i,j] * (P[i,j] - P[i-1,j])/dt +
                # 壁面换热
                q_wv/(rho_v[i,j] * alpha_v[i,j]) +
                # 相变热量
                gamma * h_v[i,j]/(rho_v[i,j] * alpha_v[i,j]) +
                # 添加燃料热生成和径向传导（气相部分）
                (q_fuel * alpha_v[i,j] + q_rad)/(rho_v[i,j] * alpha_v[i,j])
            )
            
            # ====== 12. 时间推进 ======
            new_alpha_l[i,j] = alpha_l[i,j] + dt * d_alpha_l_dt
            new_alpha_v[i,j] = alpha_v[i,j] + dt * d_alpha_v_dt
            new_u_l[i,j] = u_l[i,j] + dt * d_ul_dt
            new_u_v[i,j] = u_v[i,j] + dt * d_uv_dt
            new_v_l[i,j] = v_l[i,j] + dt * d_vl_dt
            new_v_v[i,j] = v_v[i,j] + dt * d_vv_dt
            new_h_l[i,j] = h_l[i,j] + dt * d_hl_dt
            new_h_v[i,j] = h_v[i,j] + dt * d_hv_dt
            
            # 限制物理范围
            new_alpha_l[i,j] = min(max(new_alpha_l[i,j], 0.0), 1.0)
            new_alpha_v[i,j] = 1.0 - new_alpha_l[i,j]
            
    # 在返回结果之前检查数值稳定性
    dt_new, needs_restart = check_numerical_stability(
        new_alpha_l, new_alpha_v, new_u_l, new_u_v, new_v_l, new_v_v,
        P, T_l, T_v, rho_l, rho_v, h_l, h_v, dt
    )
    
    if needs_restart:
        # 如果需要重新计算，返回特殊标记
        return (None, None, None, None, None, None, None, dt_new, True)
    
    return (new_alpha_l, new_alpha_v, new_u_l, new_u_v, 
            new_v_l, new_v_v, new_h_l, new_h_v, dt_new, False)

@njit
def calculate_turbulent_mixing_numba(alpha_i, alpha_j, rho, dy):
    """计算湍流混合项"""
    C_e = 0.01  # 涡流扩散率
    L_tur = dy  # 湍流混合长度
    alpha_eq_i = 0.5  # 平衡空泡分数
    
    return C_e * dy * L_tur * rho * (alpha_j - alpha_i - (alpha_j - alpha_eq_i))

@njit
def calculate_interfacial_force_numba(alpha_l, rho_l, u_l, u_v, v_l, v_v):
    """计算界面力"""
    C_i = 0.44  # 界面阻力系数
    d_b = 0.001  # 气泡直径
    
    # 计算相对速度
    u_r = u_v - u_l
    v_r = v_v - v_l
    V_r = np.sqrt(u_r**2 + v_r**2)
    
    # 计算界面面积密度
    a_i = 6 * (1-alpha_l) / d_b
    
    return 0.5 * C_i * rho_l * V_r * V_r * a_i

@njit
def calculate_wall_heat_flux_numba(T_w, T_f, h):
    """计算壁面热流密度"""
    return h * (T_w - T_f)

@njit
def calculate_saturation_temperature_numba(P):
    """计算饱和温度"""
    return 373.15 + 0.0001 * (P - 1.013e5)  # 简单的饱和温度计算

@njit
def calculate_latent_heat_numba(P):
    """计算汽化潜热"""
    return 2.257e6  # 常压下水汽化潜热

@njit
def calculate_cfl_timestep(u_l, u_v, v_l, v_v, dx, dy, dz):
    """计算基于CFL条件的时步长"""
    # 计算最大速度
    max_velocity = max(
        np.max(np.abs(u_l)),
        np.max(np.abs(u_v)),
        np.max(np.abs(v_l)),
        np.max(np.abs(v_v))
    )
    
    # CFL条件，取最小网格间距
    min_dx = min(dx, dy, dz)
    
    # CFL数设为0.5以保证稳定性
    cfl = 0.5
    
    # 计算时间步长
    dt = cfl * min_dx / (max_velocity + 1e-10)
    
    return dt

@njit
def check_numerical_stability(alpha_l, alpha_v, u_l, u_v, v_l, v_v, P, T_l, T_v,
                            rho_l, rho_v, h_l, h_v, dt_old):
    """检查数值稳定性并调整时间步长"""
    # 1. 检查物理量的变化率
    max_alpha_change = max(
        np.max(np.abs(np.diff(alpha_l))),
        np.max(np.abs(np.diff(alpha_v)))
    )
    
    max_velocity_change = max(
        np.max(np.abs(np.diff(u_l))),
        np.max(np.abs(np.diff(u_v))),
        np.max(np.abs(np.diff(v_l))),
        np.max(np.abs(np.diff(v_v)))
    )
    
    max_pressure_change = np.max(np.abs(np.diff(P))) / np.mean(P)
    
    max_temp_change = max(
        np.max(np.abs(np.diff(T_l))) / np.mean(T_l),
        np.max(np.abs(np.diff(T_v))) / np.mean(T_v)
    )
    
    # 2. 设置阈值
    alpha_threshold = 0.1
    velocity_threshold = 0.2
    pressure_threshold = 0.1
    temp_threshold = 0.1
    
    # 3. 根据变化率调整时间步长
    dt_factor = 1.0
    
    if max_alpha_change > alpha_threshold:
        dt_factor = min(dt_factor, 0.5)
    if max_velocity_change > velocity_threshold:
        dt_factor = min(dt_factor, 0.5)
    if max_pressure_change > pressure_threshold:
        dt_factor = min(dt_factor, 0.5)
    if max_temp_change > temp_threshold:
        dt_factor = min(dt_factor, 0.5)
        
    # 如果所有变化率都很小，可以适当增加时间步长
    if (max_alpha_change < alpha_threshold/2 and
        max_velocity_change < velocity_threshold/2 and
        max_pressure_change < pressure_threshold/2 and
        max_temp_change < temp_threshold/2):
        dt_factor = min(1.2, 1.0/dt_factor)
        
    # 4. 检查物理约束
    if (np.any(alpha_l < 0) or np.any(alpha_l > 1) or
        np.any(P < 0) or
        np.any(T_l < 273) or np.any(T_v < 273) or
        np.any(rho_l < 0) or np.any(rho_v < 0)):
        dt_factor = 0.1  # 严重违反物理约束，大幅减小时间步长
        
    # 5. 限制时间步长的变化范围
    dt_new = dt_old * dt_factor
    dt_new = max(1e-6, min(dt_new, 1e-2))  # 限制在合理范围内
    
    return dt_new, dt_factor < 1.0

@njit
def calculate_fuel_heat_generation_numba(i, n_axial, dz):
    """计算燃料热生成 - Numba优化版本"""
    q_max = 1e6  # 最大线功率密度 (W/m)
    z = i * dz
    H = n_axial * dz
    # 使用余弦分布模拟轴向功率分布
    q = q_max * np.cos(np.pi * z / H)
    return q

@njit
def calculate_radial_heat_conduction_numba(T_c, T_f, dr):
    """计算径向热传导 - Numba优化版本"""
    k = 30.0  # 假设的导热系数 (W/m·K)
    # 简化的径向热传导计算
    q_rad = k * (T_c - T_f) / dr
    return q_rad

def set_boundary_conditions(self):
        """
        设置动态边界条件
        
        包括:
        1. 入口条件的动态调整
        2. 出口条件的处理
        3. 壁面条件的设定
        """
        try:
            # 1. 入口边界条件
            for j in range(self.n_channels):
                channel_type = self.channel_types[j+1]
                
                # 1.1 计算入口空泡率
                if self.current_time < 0.1:  # 启动阶段
                    self.void_fraction[0,j] = 0.0
                else:
                    # 根据上游条件动态调整
                    self.void_fraction[0,j] = min(0.1, self.void_fraction[1,j])
                
                # 1.2 入口速度分布
                props = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')
                if props is not None:
                    # 考虑子通道类型的流量分配
                    if channel_type == 'corner':
                        factor = 0.7
                    elif channel_type == 'edge':
                        factor = 0.9
                    else:  # center
                        factor = 1.0
                        
                    self.liquid_velocity[0,j] = factor * self.inlet_mass_flux / props['density']
                    self.vapor_velocity[0,j] = 1.2 * self.liquid_velocity[0,j]  # 气相速度略大
                    
                # 1.3 入口温度和压力
                self.temperature[0,j] = self.inlet_temp
                self.pressure[0,j] = self.inlet_pressure
                
            # 2. 出口边界条件
            for j in range(self.n_channels):
                # 2.1 零梯度条件
                self.void_fraction[-1,j] = self.void_fraction[-2,j]
                self.liquid_velocity[-1,j] = self.liquid_velocity[-2,j]
                self.vapor_velocity[-1,j] = self.vapor_velocity[-2,j]
                
                # 2.2 出口压力指定
                self.pressure[-1,j] = self.outlet_pressure
                
                # 2.3 温度外推
                self.temperature[-1,j] = 2*self.temperature[-2,j] - self.temperature[-3,j]
                
            # 3. 壁面边界条件
            for i in range(self.n_axial):
                for j in range(self.n_channels):
                    channel_type = self.channel_types[j+1]
                    
                    # 3.1 壁面无滑移条件
                    if channel_type == 'corner':
                        self.liquid_velocity[i,j] *= 0.9  # 考虑强壁面效应
                        self.vapor_velocity[i,j] *= 0.9
                        
                    # 3.2 壁面温度条件
                    if hasattr(self, 'wall_temperature'):
                        self.T_c[i,j] = self.wall_temperature
                    else:
                        # 使用简化的壁温计算
                        self.T_c[i,j] = self.temperature[i,j] + 30  # 假设温差30K
                        
        except Exception as e:
            print(f"设置边界条件时出错: {str(e)}")

def apply_numerical_stabilization(self):
        """
        应用数值稳定化技术
        
        包括:
        1. 自适应时间步长控制
        2. 数值通量限制器
        3. 压力-速度耦合处理
        4. 数值振荡抑制
        """
        try:
            # 1. 自适应时间步长控制
            dt_cfl = self.calculate_cfl_timestep()
            dt_diff = self.calculate_diffusion_timestep()
            self.dt = min(dt_cfl, dt_diff)
            
            # 2. 应用通量限制器
            self.apply_flux_limiter()
            
            # 3. 压力-速度耦合处理
            self.apply_pressure_velocity_coupling()
            
            # 4. 数值振荡抑制
            self.apply_oscillation_damping()
            
        except Exception as e:
            print(f"应用数值稳定化时出错: {str(e)}")
            
def calculate_cfl_timestep(self):
        """计算基于CFL条件的时间步长"""
        try:
            # 计算最大速度
            max_velocity = max(
                np.max(np.abs(self.liquid_velocity)),
                 np.max(np.abs(self.vapor_velocity))
            )
            
            
            # 考虑声速的影响
            props_l = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')
            if props_l is not None:
                sound_speed = np.sqrt(1.4 * self.inlet_pressure / props_l['density'])
                max_velocity = max(max_velocity, sound_speed)
                
            # 计算CFL时间步长
            dt_cfl = 0.5 * self.dz / max_velocity  # CFL数取0.5
            
            return dt_cfl
            
        except Exception as e:
            print(f"计算CFL时间步长时出错: {str(e)}")
            return 1e-4  # 返回一个安全的默认值
            
def calculate_diffusion_timestep(self):
        """计算扩散稳定性时间步长"""
        try:
            # 获取流体物性
            props = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')
            if props is None:
                return 1e-4
                
            # 计算热扩散系数
            alpha = props.get('thermal_conductivity', 0.6) / \
                    (props['density'] * props.get('specific_heat', 4200))
                    
            # 计算动力粘性扩散系数
            nu = props['viscosity'] / props['density']
            
            # 取最扩散系数
            max_diff = max(alpha, nu)
            
            # 计算扩散时间步长
            dt_diff = 0.25 * self.dz**2 / max_diff  # 因子0.25确保稳定性
            
            return dt_diff
            
        except Exception as e:
            print(f"计算扩散时间步长时出错: {str(e)}")
            return 1e-4
            
def apply_flux_limiter(self):
        """应用通量限制器以抑制数值振荡"""
        try:
            for i in range(1, self.n_axial-1):
                for j in range(self.n_channels):
                    # 计算空泡率梯度
                    r_void = (self.void_fraction[i,j] - self.void_fraction[i-1,j]) / \
                            (self.void_fraction[i+1,j] - self.void_fraction[i,j] + 1e-10)
                    
                    # 使用van Leer限制器
                    phi = (r_void + abs(r_void)) / (1 + abs(r_void))
                    
                    # 限制空泡率变化
                    self.void_fraction[i,j] = self.void_fraction[i-1,j] + \
                                            phi * (self.void_fraction[i+1,j] - self.void_fraction[i,j])
                                            
                    # 对速度场应用类似处理
                    r_vl = (self.liquid_velocity[i,j] - self.liquid_velocity[i-1,j]) / \
                           (self.liquid_velocity[i+1,j] - self.liquid_velocity[i,j] + 1e-10)
                    phi_vl = (r_vl + abs(r_vl)) / (1 + abs(r_vl))
                    self.liquid_velocity[i,j] = self.liquid_velocity[i-1,j] + \
                                              phi_vl * (self.liquid_velocity[i+1,j] - self.liquid_velocity[i,j])
                                              
        except Exception as e:
            print(f"应用通量限制器时出错: {str(e)}")
            
def apply_pressure_velocity_coupling(self):
        """应用SIMPLE算法处理压力-速度耦合"""
        try:
            # 1. 求解动量预测方程
            u_star = np.copy(self.liquid_velocity)
            v_star = np.copy(self.vapor_velocity)
            
            # 2. 求解压力校正方程
            p_correction = self.solve_pressure_correction()
            
            # 3. 校正压力场
            self.pressure += 0.8 * p_correction  # 松弛因子0.8
            
            # 4. 校正速度场
            for i in range(1, self.n_axial-1):
                for j in range(self.n_channels):
                    # 速度校正
                    dp_dx = (p_correction[i+1,j] - p_correction[i-1,j]) / (2*self.dz)
                    props = self.get_fluid_properties(self.pressure[i,j], 
                                                    self.temperature[i,j], 'liquid')
                    if props is not None:
                        self.liquid_velocity[i,j] = u_star[i,j] - \
                                                  self.dt * dp_dx / props['density']
                        self.vapor_velocity[i,j] = v_star[i,j] - \
                                                 self.dt * dp_dx / (props['density'] * 0.1)
                                                 
        except Exception as e:
            print(f"应用压力-速度耦合时出错: {str(e)}")
            
def apply_oscillation_damping(self):
        """应用数值振荡抑制技术"""
        try:
            # 1. 人工粘性项
            epsilon = 0.1  # 人工粘性系数
            
            for i in range(1, self.n_axial-1):
                for j in range(self.n_channels):
                    # 对空泡率应用人工粘性
                    d2alpha = (self.void_fraction[i+1,j] - 2*self.void_fraction[i,j] + 
                              self.void_fraction[i-1,j])
                    self.void_fraction[i,j] += epsilon * d2alpha
                    
                    # 对速度场应用人工粘性
                    d2vl = (self.liquid_velocity[i+1,j] - 2*self.liquid_velocity[i,j] + 
                            self.liquid_velocity[i-1,j])
                    self.liquid_velocity[i,j] += epsilon * d2vl
                    
                    d2vv = (self.vapor_velocity[i+1,j] - 2*self.vapor_velocity[i,j] + 
                            self.vapor_velocity[i-1,j])
                    self.vapor_velocity[i,j] += epsilon * d2vv
                    
            # 2. 限制物理范围
            self.void_fraction = np.clip(self.void_fraction, 0, 0.9999)
            self.liquid_velocity = np.maximum(self.liquid_velocity, 0)
            self.vapor_velocity = np.maximum(self.vapor_velocity, 0)
            self.pressure = np.maximum(self.pressure, 1e5)  # 最小压力1bar
            self.temperature = np.clip(self.temperature, 273.15, 647)  # 水的物理范围
            
        except Exception as e:
            print(f"应用振荡抑制时出错: {str(e)}")

def solve_pressure_correction(self):
        """求解压力校正方程"""
        try:
            n = self.n_axial * self.n_channels
            A = lil_matrix((n, n))
            b = np.zeros(n)
            
            for i in range(1, self.n_axial-1):
                for j in range(self.n_channels):
                    idx = i * self.n_channels + j
                    
                    # 获取流体物性
                    props = self.get_fluid_properties(self.pressure[i,j], 
                                                    self.temperature[i,j], 'liquid')
                    if props is None:
                        continue
                        
                    # 构建系数矩阵
                    A[idx, idx] = -2
                    A[idx, idx-self.n_channels] = 1
                    A[idx, idx+self.n_channels] = 1
                    
                    # 构建源项（质量不平衡）
                    dm_in = props['density'] * self.liquid_velocity[i-1,j]
                    dm_out = props['density'] * self.liquid_velocity[i+1,j]
                    b[idx] = dm_in - dm_out
                    
            # 求解线性方程组
            p_correction = spsolve(A.tocsr(), b)
            
            # 重塑为二维数组
            return p_correction.reshape((self.n_axial, self.n_channels))
            
        except Exception as e:
            print(f"求解压力校正方程时出错: {str(e)}")
            return np.zeros((self.n_axial, self.n_channels))

def calculate_heat_transfer(self, i, j):
        """
        计算综合传热过程
        
        包括:
        1. 轴向热传导
        2. 径向热传导
        3. 壁面传热
        4. 相间传热
        """
        try:
            # 获取局部条件
            P = self.pressure[i,j]
            T_l = self.T_l[i,j]
            T_v = self.T_v[i,j]
            T_w = self.T_c[i,j]
            
            # 1. 轴向热传导
            q_axial = self.calculate_axial_heat_conduction(i, j)
            
            # 2. 径向热传导
            q_radial = self.calculate_radial_heat_conduction(i, j)
            
            # 3. 壁面传热
            q_wall = self.calculate_wall_heat_transfer(i, j)
            
            # 4. 相间传热
            q_interface = self.calculate_interfacial_heat_transfer(i, j)
            
            return q_axial, q_radial, q_wall, q_interface
            
        except Exception as e:
            print(f"计算综合传热时出错: {str(e)}")
            return 0.0, 0.0, 0.0, 0.0

def calculate_axial_heat_conduction(self, i, j):
        """计算轴向热传导"""
        try:
            if i == 0 or i == self.n_axial-1:
                return 0.0
                
            # 获取流体物性
            props_l = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
            if props_l is None:
                return 0.0
                
            # 计算温度梯度
            dT_dz = (self.T_l[i+1,j] - self.T_l[i-1,j]) / (2*self.dz)
            d2T_dz2 = (self.T_l[i+1,j] - 2*self.T_l[i,j] + self.T_l[i-1,j]) / (self.dz**2)
            
            # 计算轴向热流密度
            q_axial = -props_l['thermal_conductivity'] * d2T_dz2
            
            return q_axial
            
        except Exception as e:
            print(f"计算轴向热传导时出错: {str(e)}")
            return 0.0

def calculate_radial_heat_conduction(self, i, j):
        """计算径向热传导"""
        try:
            # 获取子通道类型和几何参数
            channel_type = self.channel_types[j+1]
            flow_area = self.channel_geometry[channel_type]['flow_area']
            wetted_perimeter = self.channel_geometry[channel_type]['wetted_perimeter']
            
            # 计算特征长度
            L_char = 4 * flow_area / wetted_perimeter
            
            # 获取流体物性
            props_l = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
            if props_l is None:
                return 0.0
                
            # 计算径向温度梯度
            dT_dr = (self.T_c[i,j] - self.T_l[i,j]) / (L_char/2)
            
            # 计算径向热流密度
            q_radial = -props_l['thermal_conductivity'] * dT_dr
            
            return q_radial
            
        except Exception as e:
            print(f"计算径向热传导时出错: {str(e)}")
            return 0.0

def calculate_wall_heat_transfer(self, i, j):
        """计算壁面传热"""
        try:
            # 获取局部条件
            T_w = self.T_c[i,j]
            T_l = self.T_l[i,j]
            T_v = self.T_v[i,j]
            P = self.pressure[i,j]
            
            # 计算饱和温度
            T_sat = PropsSI('T', 'P', P, 'Q', 0, 'Water')
            
            # 根据传热机制选择合适的换热系数
            if T_l < T_sat and T_w > T_sat:
                # 过冷沸腾
                h = self.calculate_subcooled_boiling_htc(T_w, T_sat, T_l, P)
            elif T_l >= T_sat and T_w > T_sat:
                # 饱和核态沸腾
                h = self.calculate_nucleate_boiling_htc(T_w, T_sat, P)
            else:
                # 单相对流换热
                h = self.calculate_single_phase_htc(T_w, P, 'liquid')
                
            # 计算壁面热流密度
            q_wall = h * (T_w - T_l)
            
            return q_wall
            
        except Exception as e:
            print(f"计算壁面传热时出错: {str(e)}")
            return 0.0

def calculate_interfacial_heat_transfer(self, i, j):
        """计算相间传热"""
        try:
            # 获取局部条件
            T_l = self.T_l[i,j]
            T_v = self.T_v[i,j]
            P = self.pressure[i,j]
            alpha = self.void_fraction[i,j]
            
            # 获取流体物性
            props_l = self.get_fluid_properties(P, T_l, 'liquid')
            props_v = self.get_fluid_properties(P, T_v, 'vapor')
            
            if props_l is None or props_v is None:
                return 0.0
                
            # 计算界面面积密度
            d_bubble = 0.001  # 假设气泡直径为1mm
            a_i = 6 * alpha / d_bubble
            
            # 计算相对速度
            v_r = abs(self.vapor_velocity[i,j] - self.liquid_velocity[i,j])
            
            # 计算界面换热系数
            Re_i = props_l['density'] * v_r * d_bubble / props_l['viscosity']
            Pr_l = props_l['viscosity'] * props_l['specific_heat'] / props_l['thermal_conductivity']
            Nu_i = 2 + 0.6 * Re_i**0.5 * Pr_l**0.33
            h_i = Nu_i * props_l['thermal_conductivity'] / d_bubble
            
            # 计算相间传热量
            q_interface = h_i * a_i * (T_v - T_l)
            
            return q_interface
            
        except Exception as e:
            print(f"计算相间传热时出错: {str(e)}")
            return 0.0

def update_temperature_fields(self):
        """更新温度场"""
        try:
            for i in range(1, self.n_axial-1):
                for j in range(self.n_channels):
                    # 获取各传热分量
                    q_axial, q_radial, q_wall, q_interface = self.calculate_heat_transfer(i, j)
                    
                    # 获取流体物性
                    props_l = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
                    props_v = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')
                    
                    if props_l is None or props_v is None:
                        continue
                        
                    # 计算温度变化率
                    dTl_dt = (q_axial + q_radial + q_wall + q_interface) / \
                             (props_l['density'] * props_l['specific_heat'])
                         
                    dTv_dt = -q_interface / (props_v['density'] * props_v['specific_heat'])
                    
                    # 更新温度
                    self.T_l[i,j] += self.dt * dTl_dt
                    self.T_v[i,j] += self.dt * dTv_dt
                    
                    # 更新壁温
                    self.update_wall_temperature(i, j)
                    
        except Exception as e:
            print(f"更新温度场时出错: {str(e)}")

def update_wall_temperature(self, i, j):
        """更新壁面温度"""
        try:
            # 获取子通道类型和几何参数
            channel_type = self.channel_types[j+1]
            heated_perimeter = self.channel_geometry[channel_type]['heated_perimeter']
            
            # 计算燃料产热
            q_fuel = self.calculate_fuel_heat_generation(i, j)
            
            # 计算壁面导热
            k_wall = 16.0  # 假设壁面材料为不锈钢
            delta_wall = 0.001  # 假设壁厚1mm
            
            # 计算新的壁温
            self.T_c[i,j] = self.T_l[i,j] + \
                            (q_fuel * heated_perimeter * delta_wall) / k_wall
                        
        except Exception as e:
            print(f"更新壁面温度时出错: {str(e)}")

def calculate_fuel_heat_generation(self, i, j):
        """计算燃料产热"""
        try:
            # 使用余弦分布模拟轴向功率分布
            z = i * self.dz
            q_max = 1e6  # 最大线功��密度
            
            # 考虑子通道类型的功率因子
            channel_type = self.channel_types[j+1]
            if channel_type == 'corner':
                power_factor = 0.8
            elif channel_type == 'edge':
                power_factor = 0.9
            else:  # center
                power_factor = 1.0
                
            # 计算局部功率密度
            q = q_max * power_factor * np.cos(np.pi * z / self.channel_length)
            
            return q
            
        except Exception as e:
            print(f"计算燃料产热时出错: {str(e)}")
            return 0.0

@njit(parallel=True)
def compute_grid_points_numba(n_axial, n_channels,
                            void_fraction, prev_void_fraction,
                            vapor_velocity, liquid_velocity,
                            prev_liquid_velocity, prev_vapor_velocity,
                            pressure, temperature,
                            time_step, dz):
    """使用Numba优化的网格点计算"""
    new_void = np.zeros_like(void_fraction)
    new_vl = np.zeros_like(liquid_velocity)
    new_vv = np.zeros_like(vapor_velocity)
    
    for i in prange(1, n_axial-1):
        for j in prange(n_channels):
            # 计算空间导数
            d_alpha_dx = (void_fraction[i+1,j] - void_fraction[i-1,j]) / (2*dz)
            d_vl_dx = (liquid_velocity[i+1,j] - liquid_velocity[i-1,j]) / (2*dz)
            d_vv_dx = (vapor_velocity[i+1,j] - vapor_velocity[i-1,j]) / (2*dz)
            
            # 计算压力梯度
            dp_dx = (pressure[i+1,j] - pressure[i-1,j]) / (2*dz)
            
            # 更新空泡率
            new_void[i,j] = prev_void_fraction[i,j] - time_step * (
                vapor_velocity[i,j] * d_alpha_dx
            )
            
            # 更新液相速度
            new_vl[i,j] = prev_liquid_velocity[i,j] - time_step * (
                liquid_velocity[i,j] * d_vl_dx +
                dp_dx / 1000.0 +  # 假设液相密度约为1000 kg/m³
                9.81  # 重力项
            )
            
            # 更新气相速度
            new_vv[i,j] = prev_vapor_velocity[i,j] - time_step * (
                vapor_velocity[i,j] * d_vv_dx +
                dp_dx / 1.0 +  # 假设气相密度约为1 kg/m³
                9.81  # 重力项
            )
            
            # 限制物理范围
            new_void[i,j] = min(max(new_void[i,j], 0.0), 0.9999)
            new_vl[i,j] = max(new_vl[i,j], 0.0)
            new_vv[i,j] = max(new_vv[i,j], 0.0)
            
    # 处理边界条件
    new_void[0,:] = void_fraction[0,:]
    new_void[-1,:] = void_fraction[-2,:]
    new_vl[0,:] = liquid_velocity[0,:]
    new_vl[-1,:] = liquid_velocity[-2,:]
    new_vv[0,:] = vapor_velocity[0,:]
    new_vv[-1,:] = vapor_velocity[-2,:]
    
    return new_void, new_vl, new_vv

if __name__ == "__main__":
    # 创建分析实例
    analysis = SubchannelAnalysis(n_axial=20)
    
    # 运行瞬态模拟
    analysis.run_transient_simulation()
    
    # 显示瞬态结果
    analysis.plot_transient_results()
    
    # 显示最终状态的详细分布
    analysis.plot_results()